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Abstract 


The method of weighted residuals in the form ">f a modified Galerkin 
method with boundary residuals is developed for the study of the transmission 
of sound in nonuniform ducts carrying a steady, compressible flow. In this 
development the steady flow is modeled as essentially one dimensional but with 
a kinematic modification to force tangency of the flow at the duct walls. 

Three forms of the computational scheme are developed using for basis functions 
(a) the no-flow uniform duct modes, (b) positive running uniform duct modes, 
with flow, and (c) positive and negative running uniform duct modes, with 
flow. The formulation using the no-flow modes is the most highly developed and 
has advantages primarily due to relative computational simplicity. Results 
using the three methods are shown to converge to known solutions for several 
special cases. The most significant check case is against low frequency, one 
dimensional results over the complete subsonic Mach number range. Development 
of the method is continuing, with emphasis on assessing the relative accuracy 
and efficiency of the three implementations. 
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INTRODUCTION 


The method of weighted residuals (MWR) in the form of a Galerkin 
technique with boundary residuals, has been shown to be a useful method 
for the investigation of the transmission of sound in nonuniform ducts 
with no flow^^. In this paper the method is extended to treat the 

case when the duct carries a steady nonuniform flow. 

The introduction of steady mean flow into the study of nonuniform 
duct propagation complicates the problem in two ways. The acoustic 
field equations no longer reduce to either the simple wave equation of 
the noflow case or the convective wave equation of the uniform flow 
case. In addition, the field equations represent acoustic perturbat- 
ions on a steady nonuniform flow field which in itself is difficult to 
describe. As a consequence of these complications, no simple theories 
of propagation have been developed and investigations appearing in the 
literature have been relatively few. 

The majority of studies appearing have been based on a one dimension- 
al theory which treats the propagation as plane waves moving in a one 

( 2 ) ( 3 ) 

dimensional nozzle flow. Powell , Eisenberg and Kao King and 
Karamcheti ^ , Huerre and Karamcheti ^ and Davis ^ have studied one 
dimensional models by a variety of techniques. 

It is the purpose of this investigation to study multi-modal prop- 
agation in nonuniform ducts with flow. Published work dealing with 
this problem has been limited to approximate methods for nonuniformities 

in cross section, lining and flow properties which are slowly varying. 

( 7 ) 

Tam considered hardwalled ducts with slowly varying cross section by 

using a Born type of approximation and Fourier Transform methods. Nayfeh 

( 8 ) 

and. his co-workers have used the method of multiple scales and include 
slowly varying area, lining and boundary layer. While Tam’s work 

allows modal coupling, Nayfeh’ s does not. 

The method described here treats the complete acoustic field equations 
and is approximate only in the sense that the solutions are represented 
in terms of a superposition of a limited number of specified functions. 

We will be concerned primarily with the method of solving the field 
equations and relatively little with the difficulty presented by describ- 
ing the nonuniform steady flow field. Instead, we will use an approx- 
imation for the flow field based on one dimensional compressible flow 
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with a kinematic modification to allow for flow tangency at the duct walls. 

In addition, in this initial study of computational schemes we will 
concentrate primarily on the solution in elementary duct sections. Rather than 
an exhaustive treatment of the transmission characteristics of various duct 
configurations we limit our attention to the successful implementation of the 
scheme in particular cases. 

A principal motivation for the study of propagation in ducts with flow 
is the observed phenomenon of subsonic acoustic choking which occurs in inlet 
type flows when the acoustic source radiates upstream through a flow constriction. 
It is found that there is a substantial reduction in acoustic transmission 
past the constriction when the Mach number in the constriction exceeds about 
M = 0.7. One would predict this intuitively if the constriction is sonic; 
however, the occurence of acoustic choking at such low Mach number is surprising 

(9) 

and requires investigation. Tam has considered the possibility that some 
of the choking is attributable to nonlinear effects arising when high intensity 
sound propagates upstream. The computational scheme developed in this paper 
provides a basis for studying the linear aspects of this problem. Extensive 
computations directed toward this goal will be made in an extension of the 
current research effort. 

CO-ORDINATE SYSTEM AND GEOMETRY 


In the following development we will consider only two dimensional ducts. 
The method employed generalizes to the three dimensional case, but is then, 
of course, more computationally demanding. We will also consider e^qglicitly 
only ducts of infinite length, although as has been shown ^ , the extension 
to ducts which are of finite length is not difficult provided that the reflection 
properties of the termination are known. 

Figure 1 shows the type of configuration under consideration in which 
two semi-infinite uniform duct sections of height b^ and b^ are joined by a 
transition section of variable height b(x). The duct wall at y = 0 is hard and 
the wall at y - b(x) in the transition section has variable impedance Z o (xj . 

The uniform sections have impedance and Z^ at y = b Q ana y = b^, respectively, 

and have a hard wall at y = 0. 
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This representation can be considered as a model of a duct of this con- 
figuration or of one symmetric with respect to the x axis with symmetric 
propagation. We will consider only continuous variations in duct 
height; however discontinuous changes in impedance will be permitted at 
the ends of the nonuniformity. 
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GOVERNING ACOUSTIC EQUATIONS 


We will consider acoustic propagation as small perturbations on the 
steady duct flow. It is assumed in the present analysis that the fluid 
motion is nonviscous and isentropic. To derive the acoustic equations 
we begin with the equations of continuity, momentum, energy and state 
in dimensional form 

-h dnr f'r* = ° 
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where p* , p* , v* are the steady flow conditions defined by 


dir 

• O 

(e) 

\£ — — JL£ 4 *. 

_ jjKad J?c, ¥ 

(3) 

V o *. Jr a ef 

+ at ns t/ 0 ¥ - o 

(to) 

fa 

¥ yX 1 

(») 

fr 

■/ 


Because of the isentropic assumption the energy equation is directly deriv- 
able from the continuity equation * The energy equation does not contain 
the density and is more convenient to use. By using Equations (5) - (7) 
in Equations (2) - (4) and by making the acoustic assumption of small 

perturbations on the steady flow, the dimensional acoustic momentum, energy 
and state equations are: 
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Equations (9) and (10) were used to eliminate the steady flow terms in 
Equations (12) and (13) . The term in Equation (12) containing p* can be 
rewritten by using Equation (14) and equation (9) : 
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The modified momentum equation is then 
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The governing acoustic equations can be non-dimensionalized in the 
standard way by defining the nondimensional variables. 


7°~~ 




/*" 


/*<- 


/* 


c- 






* 


O- 


The reference states^ , c r are arbitrary values of the density and speed 
of sound. They will generally be chosen as the state which exists in the 
uniform flow incident on the nonuniformity. The nondimensional equations 
of momentum, energy and state are then 
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is assumed, Equations 
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When a harmonic time dependence of the form e 
(16) and (17) become 
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where Jc r = ao/c r - nr/A.,. and k, r i s the free space wave number 

in the reference state. 


In order to specify a boundary condition at the duct walis, we introduce 
a co-ordinate system as shown in Figure 1. This figure shows the manner 
in which the duct height profile and slope of the height profile are 
specified as well as the designation of the local outward unit normal 
at the duct wall, y\ The duct wall boundary condition employed is charac- 
teristic of a normally reacting lining in the presence of a harmonic 
pressure variation. 






where v*^ is the particle velocity in the acoustic lining at the surface , 
assumed in the direction of the outward normal. Z is the wall impedance 
which may be a function of axial position. In nondime nsional variables the 
boundary condition becomes 
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or 
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where A = p c /2 is the acoustic admittance ratio of the lining based on the 
r r 

reference admittance 1/p c . In harmonic motion the relation between the 

r r 

component of fluid particle velocity at the wall and particle velocity in 
the lining at the wall is given by 

v r -y‘ 
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U is the liurd velocity tangent to the wall and 
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derivative along the wall. In nondimensional form this becomes 
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The boundary condition at a duct wall, assuming a time dependence e , 
is thus 

7-y* Jr- 1 jfr ( A 7=) 

Equations (19) and (20) with the duct wall boundary condition of 
Equation (21) specify the boundary value problem for propagation in a 
variable geometry duct with a steady flow. The problem as described is 
formidable, principally due to the introduction of the steady duct flow. 

The introduction of the flow field considerably complicates the acoustic 
problem, but it also presents the requirement of describing the flow field. 
The flow field would be defined by the solution of Equations (8) - (11) 
pi ^3 suitable boundary conditions. Even for simple duct geometries this 
requires sophisticated numerical techniques and only rarely in highly 
specialised cases co-id one hope to generate an exact closed form solution. 
For detailed studies a fairly exact description of the flow field may be 
required. However, in order to study the effect of duct nonuni formi ties 
on the propagation of sound with a certain degree of generality, we will 
use an approximation to the steady flow field. 


The most common approximate description of the steady flow field is 
the one-dimensional theory which is one of the cornerstones of elementary 
gasdynamics. This theory proves very useful for a wide range of duct 
contours, but with rigor must be viewed as a first approximation when the 
area change is gradual. 


In the one-dimensional theory the variation of the Mach number, 
nondimensional density and pressure with axial position are given by 
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where A is the local cross sectional area, 
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In addition, the local nondimensional speed of sound is given by: 
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and the local nondimensional flow velocity is given by: 
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U q is the local axial component of flow, the transverse components being 
assumed zero in the one dirensional approximation. 


In this approximation the flow properties are assumed to be uniform 
at a cross section with no transverse velocity. They are solutions of the 
one dimensional continuity and momentum equations 
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together with the isentropic equation of state 



The acoustic problem. Equations (19) and (20 ) , involves axial and 
transverse velocities and pressures as well as axial and transverse gradients 
Furthermore, the boundary condition requires the velocity tangent to the wall 
This required velocity field could be approximated simply by using the one 
dimensional theory and ignoring the transverse velocity and the transverse 
gradients. We choose to include some measure of the effect of transverse 
velocity and its gradient by taking the solution to Equations (8) - (11) to 
be, in nondimensional form for the two dimensional duct 

Uo (x, {/) = U 0 fx) 
v.(x.cf)-- o.(x>Ji 

c a (x,y) - c 0 (x) 
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where a bar notation has been used to denote the one dimensional solution. 
This form for the steady flow field must be recognised as approximate, but 
it is considered to be a reasonable approximation in that it exactly satis- 
fies both the continuity equation. Equation (8) and the x component of the 
exact momentum equation. Equation (9) . The y component of Equation (9) is 
not satisfied exactly, but is satisfied on an average basis on the cross 
section of a duct symmetric with respect to the x axis. In addition, this 
solution satisfies the requirement of flow tangency at the duct wall. 

By using this approximation we are introducing an error into Equations 
(191 and (20) since in those equations it is assumed that the flow field 
satisfies Equations (8) - (11) . Based on the success of one dimensional 
nozzle theory, the axial velocity, pressure and density variations should 
be acceptable over a fairly wide range of duct shapes. The assumption for 
the transverse velocity is on a less firm basis , and will be more accurate 
for gradual T .rea changes. In the absence of detailed information about 
the flow licid it is felt that tha approximate flow field will be usefu] 
for identifying the important general properties of transmission in non- 
uniform ducts. 


The steady pressure, density and axial velocity in the steady flow 

approximation are functions of x alone so that the steady flow field is 

described by U (x) , V (x,y) , p (x) , p (x) . In addition the transverse 
o o o o 

gradient 3V Q /3y is a function of x alone. These simplifications provide 
significant reductions in numerical computations. However, if more exact 
flow field descriptions are used these restrictions must be relaxed. 


With this approximation for the flow field, the momentum and energy 
equations, Equations (19) and (20) can be expanded in the present two 
dimensional case to yield: 
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The boundary condition of Equation (21) can be expanded by noting that 
at the duct wall: 

U T = U a CO S 9 V- V Q s >s> & 
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The boundary conditicr: at the wall for the steady flow is that the normal 
velocity must vanish: 
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Equation (21) becomes, at a duct wall 
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METHOD OF SOLUTION IN THE NONUNIFORM SECTION 


In Reference (1), the method of weighted residuals (MWR) (modified 
Galerkin method) , was employed to study acoustic propagation in nonuniform 
ducts without flow. We will use this method with minor changes in the 
present study. (Reference (l) is reproduced in Appendix (A)). 
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Attention will be restricted to two dimensional ducts with geometry 
as shown in Pig. 1. As previously noted, this can be considered as a 
model for a duct actually shaped as shown or else one symmetric with respect 
to the x axis. In the latter case we consider only symmetric propagation. 

We seek solutions to the field equations, Equations (22) - (24), and 
the boundary condition, Equation (25) in the form: 

■pac'tj)* Z r„ m ( *' y) 

n 

U (X, <j) = Z <*> 6 » 

f? 

v y) = Z (x) </>„ ( X, y) 

The success of the MWR is dependent on an appropriate choice of basis 
functions. In the no flow case of Reference (1) , the basis functions were 
chosen to he the transverse modal functions which would exist in a 
uniform duct with the properties existing locally in the nonuniform duct. 
This philosophy can be carried over to the duct with a steady mean flow. 

In this case the appropriate choice of basis functions would be the trans- 
verse modal functions for a uniform duct with flow and geometry properties 
existing locally in the nonuniform duct. In the nonuniform duct the flow 
properties vary transversely so that the equivalent uniform flow to be used 
in defining the basis functions is open to interpretation, since the 
boundary residual will be more nearly satisified if the basis functions 
nearly satisy the boundary condition, it would seem appropriate to use 
the velocity at the wall as the equivalent uniform velocity. Of course 
in the simple flow model used here this is unimportant, since the wall vel- 
ocity is the same as the velocity on the duct axis. 

The major disadvantage to using this type of basis function is that 
there are different modes for upstream and downstream propagation so that 
to obtain the same resolution as in the no flow case we would use twice as 
many modes, half representing upstream propagation and half representing 
downstream propagation. 

A second possibility would be the use of tne no flow modes. Since 
these modes are not generally near-solutions we would expect to generate 
solutions which require a large number of terms to converge. We tried 
using the no flow modes in the initial development and found that they 
worked well at low hach numbers with soft walls and ^or nearly hard walls. 
In both instances . the., no flow and flow tr ansverse modal functions are not 
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much different so that this success is not unexpected. For this reason 
we have one version of the computational scheme for nearly hard walls or 
for soft walls at low Mach number. We will detail here the more general 
case. 


Using the basis functions from uniform duct theory {with flow) , 

solutions are sought in the form: 
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The Mach number to use is that which exists at the duct wall. This 
eigenvalue problem simultaneously yields the transverse wave number 
and k for a mode of propagation. Eigenvalues for propagation in 

X 

n 

both the plus and minus axial direction will appear and will be identified 

by the sign in k . 

n 

The eigenfunctions which arise i 1 ' Lhe flow case are not orthogonal 
and the orthogonality property cannot be used in computations. However, 
in the flow case this disadvantage is not of serious consequence in view of 
the other complexity which already exists. 

If the assumed solution is substituted in the differential equations, 
Equations - (24), it will in general not satisfy them and the left 

hand side will not equate to zero. Instead, it will add up to an error, 
or residual. By denoting the residual by R, we w’-ite upon substituting 
the assumed solution in the governing equations: 


13 . 


R, = u a + _L 


3x 


/* 


3X 


+ \/ 0 


+ 


(i k r +■ - -L fP 

K r £>*. s " **) /) / v 




& z = u 0 ^f* + y ay* j_ ay* yp u ( l ^ + ay 0 

^ <5>i/ /2 ay ?x r * 


ax 


R 3 - tflfg. , U'ir* + V.&. + tifi, + 2 £ Un 

+ y+*(lg + Tp]r« 


In a similar manner we form a boundary residual, 
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complete set. (See Eversman, Cook and Beckemeyer ^ , Finlaysan and 
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A unit depth of the duct is assumed. Note that the trial functions 
already satisfy the boundary condition at the hard wall (or duct centre 
line) so that no boundary residual is required there. In expanded form 
the orthogoniality conditions for the differential equations yield 
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Equations (26), (27), (28) can be rewritten and simplified to some 
extent by following the sequence of steps outlined in Reference (1) (See 
Appendix A of the present document) . We use rules for differentiation 
of integrals, integration by parts, and the condition of flow tangency 
at the wall to obtain new forms of Equations (26) and (27) for n = 1,2,...N: 
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The weighted boundary residual takes the form 
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By using the steps outlined above, plus the boundary residual. Equation (28) 
is rewritten for n = 1,2,...N 
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The Method of Weighted Residuals can be completed by inserting in 
Equations (29) , (30) , and (31) the truncated serias representations for 
p , u N , v . If this is done and the terms are grouped there results 
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Equation (34) required special consideration in that the terms contributed 
by the boundary residual yield derivatives of (this did not occur in the 
no-flow case) . We have introduced the definitions 
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We make use of the definitions of the basis functions 
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the d rivative relations 
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Using these relations, we arrive at definitions of the coefficients in 
Eqs (35) : 
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Equations (35) can be conveniently written in matrix form 
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In order to achieve a form more closely related to the no flow equations 

we divide the equations by YPq/ multiply the equations by p Q and 

divide the v equations by U_. 
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MATCHING OF THE NONUNIFORM SECTION TO UNIFORM SECTIONS 


As illustrated in Figure 1, the nonuniform duct segment is considered to 
be a transition section between two uniform, infinite ducts. In this section 
we establish the procedure for matching the nonuniformity to the uniform sec- 
tions . 


The set of linear ordinary differential equations for the axial variation 
of the coefficients in the assumed solution, given by Equation (36), can be 
represented in the form 


f 



The dimensioning of the matrices involved depends on the types of basis 
functions used. In our investigations to date we have used three different 
types of basis functions; 


(1) Initial studies utilized as basis functions the modes from the no- 
flow problem. In this case the assumed solution and the required 
eigenvalues are exactly as described in Reference (1) which is 
included in Appendix A to this document. The differential equation 
is then dimensioned as follows: 






«f / 


N is the number of basis functions used. Note that in this form 
the problem is of size 3N x 3N as compared to 2N x 2N in the no 
flow case. Extensive computations were carried out using this 
method and it was found to give results in generally good agree- 
ment with known correct results for low Mach numbers with soft 
walls and for a more extensive Mach number range with nearly hard 
walls. Computations using this approach form a large portion of 
the results generated to date and are discussed in a subsequent 
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section. As will be seen, certain inadequacies exist and we have 
extended the computational scheme to use modes more closely related 
to the flow problem. The details of this approach will not be 
given here as they can be inferred from the more sophisticated 
approaches to follow. 

(2) Since the success of Gaierkin type methods depends to a large extent 
on an appropriate choice of basis functions we have modified the 
basic computational scheme described above to use modes from uniform 
duct theory with the effects of flow included. In the first step 

in this direction we employed the flow modes but included only modes 
corresponding to propagation in the positive axial direction (in the 
flow case, as distinct from the no-flow case, there is one set of 
positive moving modes and another set of negative moving modes) . 

This is really a sub-case of the more general use of the flow modes 
described below and as such will not be described separately. The 
dimensioning in this case is the same as the case w^en no-flow modes 
are used, however, the eigenvalues associated with the basis 
functions require a more extensive eigenvalue routine, as described 
in Appendix B . 

(3) The most advanced version of the Gaierkin method developed to date 
uses the full set of uniform duct flow modes, including both positive 
and negative running waves. The principal reason for avoiding this 
approach at the outset was the nominal doubling of the number of 
modes required. If we let N denote the number of positive running 
modes, then, since for each positive running mode there is a 
corresponding negative running one, we actually require 2N modes to 
achieve the same level of resolution used in the previously discussed 
cases. In this instance the dimensioning is 



btJxi 
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As noted, N is the number of inodes in one direction and is a measure 
of the basis function resolution available, consistent with the other 
implementations. At the time of preparation of this document 
numerical experimentation to verify and improve this formulation 
is still in progress. Certain new computation problems have been 
observed and are currently under investigation. Hence, the use of 
the complete flow modes can at this stage only be considered in the 
development stage. However, since this formulation is the most 
general, we will use it to explain the matching procedure. The 
corresponding procedure in the simpler formulations can then be 
easily inferred. 


The solution to the matrix differential equations of Eq. {36 ) can be 

given in terms of a transfer matrix relating u , p , v at x ■ l to u , 

* n n n n 

p , v at x = 0: 
n n 


U 


M \f. 






lYt 


' n ^ J x=o 

The transfer matrix is readily obtained by a fourth order Runge-Kutta scheme. 
We have experimented with other schemes and consider it quite probable that 
in future development some advantage, particularly in speed, may be gained 
by using a different integration schema. 
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The propagation in the uniform ducts x < 0 and x > l can be expressed in 
terms of the classical duct Lheory. In a uniform duct u£x,y), p(x,y), v(x,y) 
can be written 
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The plus and minus superscripts denote right (positive) and left (negative) 

moving modes. The K are the eigenvalues in the uniform duct. The 

HI 

determination of the direction of propagation of a given mode is discussed 
in Appendix B in the general case of softwall ducts. We have generally 
carried out our computations treating the uniform ducts as hardwalled. In 
this case 


K t = ) m-' /, Z, ■■■ v 

The choice of sign in for positive and negative moving modes can then 

(13) 

be deduced by energy flow arguments 

The nondimens ion al speed of sound c= c*/c and density p = p*/p appear 

r o o r 

because the classical duct solution used here is based on nondimensionalization 

with respect to the duct density and speed of sound, c* and p*. These 

reference conditions will generally be different in the ducts x < 0 and 

x > £. In the duct x < 0. c = 1 and p - 1 because we have defined c 

o r 

and as conditions existing in the duct x < 0. In the duct x > 0 we 

define c = c p and ■ p ? , the nondimensional speed of sound and density at 
the end of the nonuniformity. 
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The velocities and pressure in the nonuniformity can be written in matrix 
form 
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The eigenvalues in the nonuniform section are ordered so that m = l,2,.. r N 
correspond to right running modes and m ~ N-t-1, N+2,...2N correspond to left 
running modes. The method of obtaining the < and the technique used to sort 
them for the direction of propagation is discussed in Appendix B in a print 
of a paper to appear in the Journal of Sound and Vibration. 


To demonstrate the matching procedure we write £q. (39) for x = 0 
and x = SL in abbreviated form 
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As described in Reference (1) , the amplitude coefficients at x 

_ t Ax J? - 1 

absorbed the exponential terms Q m , £? 

thus amplitude coefficients referenced to x = 1 , rather than x 

(40) is written in abbreviated form at x = 0 and x = 4 
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At x = O f we match u,p,v to obtain 


( 44 ) 
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We use a least squares estimate for the u , p , v . This is obtained in 
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the form 


where 


and 


(“*, } 

' > - 
/ tv 

V*. 

x= 0 

[* eo> ]- 

. 

6MX2SJ 

■* 

<1 

a~„ 

l Vx l 

6 tit/ 



M = 

[cc.] 

[ c d °l 

x-ZA/ 

iHxi/ V 

&A/X2SS 


r h ~ 

~T~ 

l*c. ] - 

/ [Co'_ 

O 

j 

£ 


O 6/^x3 

3x6a/ 

M = 

T 

/ [Co ] [ <ol d 3 

C=N X ZA/ 

J 

O 

6N*3 

3x2-^ 


(4S) 


T 


The operation [0 o J denotes the complex conjugate transpose. This 

Put-ation can be viewed as establishing the starting values for u , p 

mm 

from the wave amplitudes in the duct x < 0. 


com- 


v 


m 
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At x = H, we match u, and p: 
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The matrix [On] is truncated to delete the 2N columns related to v . 

m 

again use a least squares matching to obtain 


We 


where 
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Again, the operation L^bj] 7 signifies the complex conjugate transpose. 

We use this notation with the general case in mind when the uniform duct 
is softwalled. When it is hardwalled, as generally assumed in our 
computations, we have dispensed with the conjugate because the eigenfunctions 
are orthogonal and a straight Fourier matching is possible. 


This operation is viewed as determining the wave structure in the 

ductx > i from the conditions at the end of the nonun if ormity. We are 

considering the possibility of matching u,p,v at x = l, rather than just 

u and p. The extra computation is insignificant and might prove beneficial 

* ( ■}* - 

since the extra information available for determining c and c might 

J m m v 

help smooth out any inaccuracies which have been introduced in the process 
of developing the transfer matrix. This procedure has not as yet been 
implemented. 
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Equations (38), (45) and (46) can now be combined to yield 
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The matrix [ ] is truncated with the 2N rows corresponding to the v equations 

deleted. 


From this point on the technique of obtaining the matrices of reflection 
and transmission coefficients follows Reference (1) . We assume an infinite 
duct x > l. Then 



From this we obtain 


faj = - [ner] {a*} us) 
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where jTrEfJ and [teanJ are the reflection and transmission coefficient 
matrices, respectively. The extension to finite ducts x > l follows the 
development of Reference (1). 



31 . 


RESULTS 


As noted previously, the development of the computational scheme has gone 
through three stages, beginning with the use of no-flow inodes for basis 
functions, progressing to positive running flow inodes as basis functions, and 
then finally to the use of the complete set of flow modes as basis functions. 

Each st ace of development represents a new level of complexity. When no- flow 
modes are used the formulation follows almost directly from the no- flow case 
described in Reference (3.) with the use of a new set of governing equations, as 
described by Eqs. {36) of this document. In advancing to the use of the 
positive propagating flow inodes the eigenvalue scheme of Reference (1) must 
be abandoned for the one described in Appendix B. Both of these implementations 
represent a size increase over the equivalent no- flow problem so that 
computational time approximately doubles. The ultimate level of complexity 
arises in using both positive and negative running flow modes as basis functions. 
Assuming that the same level of resolution is available by using basis functions 
of similar character this formulation leads to computational times in the 
neighborhood of nine times that of the original no flow case. It is apparent 
that a very careful evaluation with regard to accurracy and computing cost must 
be made to ascertain the appropriate level of complexity required to treat the 
problem. 

The major problem faced in this program has been the almost total lack of 
any results against which to make comparisons and evaluation. As noted in 
Reference (1) , we were forced to develop several alternative computational 
schemes to make an evaluation of the results of the no- flow computational 
scheme. None of the alternative schemes of the no-flow case are available 
to us in the flow case. The only approach we have been able to take is the 
reduction of the general scheme to certain special cases. Of particular interest 
in this regard is the computational scheme of Davis and Johnson in which they 
treat plane wave propagation through a one dimensional compressible flow in a 
variable cross section duct. Their computational scheme should allow us to 
generate results which serve as low frequency comparisons for the computational 
schemes developed here. 


Other comparisons that we have used serve only as gross indications of the 
consistency of our present formulation. We have approached zero flow speed 
and compared the results against our original no flow results. Since the 
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governing equations in the flow case don't degenerate directly to the no flow 
equations this type of comparison provides an independent check of some of 
the terms. We have also used the straight hardwalled duct as a comparison since 
we know that reflection coefficients should vanish and that transmission 
coefficients can be computed by simple hand computations. A further simple 
test case considers uniform soft wall ducts, In this case we are again able 
to compare against simple hand calculations. Finally, we have generated 
results for a uniform soft wall segment between two infinite hardwalled 
ducts. These results have been compared against similar results generated by 
another investigation. 

In this section we will discuss our computations to date in the form of 
a series of comparisions as noted above. The mo^t extensive results so far 
obtained are for the simplest case in which the no-flow basis functions are 
used. Since the use of flow mode basis functions is a recent development, 
only a less extensive series of results is available. In particular, as we 
will discuss later, the use of both positive and negative propagating modes 
has lead to a series of new computational problems not encountered previously. 

As a consequence our ma^or effort in this case has been directed toward isolating 
the problems and taking corrective measures. 

A, COMPUTATIONAL COMPARISONS USING NO-FLOW BASIS FUNCTIONS AND POSITIVE RUNNING 
FLOW BASIS FUNCTIONS 


As we have indicated, the use of no-flow {NF) basis functions or positive 
running flow (PF) basis functions represents a large savings in computational 
complexity and cost in comparison to the use of the complete set of positive 
and negative running flow (PNF) basis functions. Numerical experimentation 
has been carried out to determine if the NF and PF basis functions provide 
satisfactory results. 


The NF basis func.tions are certainly the simplest to use since the 
associated eigenvalues follow from a simple eigenvalue equation and it seems 
fairly certain that the functions form a complete set, each member of the set 
being orthogonal to the others. The completeness is essential for the Galerkin 
approach and the orthogonality is an advantage computationally since the 
coupling between modes in minimized. One form of our computational scheme uses 
the NF modes. As a step toward using PNF basis functions we have developed 
a form of the computer program which uses the PF basis functions. This 
formulation is not much more complicated than when the NF basis functions are 
used, except that the associated eigenvalue equation is more difficult to use. 


33 . 


The basis functions in this case are not orthogonal and we observe more coupling 
between the modes. Of particular concern is the completeness of the set of 
basis functions. Since the PF basis functions are a subset of the PNF modes, 
they may not themselves be complete. At the current stage of development our 
only means of ascertaining their usefulness i n numerical experimentation. This 
experimentation has been done in the following test cases: 

(a) Convergence of results at low Mach numbers to results obtained with 
the original no- flow duct program. 

(b) Convergence to the straight hardwall duct case in which hand cal- 
culations can be made. 

(c) Convergence to the straight uniform softwalled duct case in which hand 
calculations can be made. 

(d) Comparison of results of several schemes for softwalled, straight 
segment between two hardwalled, infinite ducts, 

(e) Comparison against a one dimensional, low frequency approximation. 

The following sub-sections describe the comparisons made to date. 

(a) Convergence to no flow case 

As reported in Reference (1) , the computational scheme applicable when no 
mean flow is present has been highly successful and provides a convenient 
base-line against which to check the results of the present program. When flow is 
present the governing equations, as given by Eqs. {36) , are not computationally 
the same as the no flow case, even when the flow Mach number becomes small. To 
be made computationally the same a number of reductions and substitutions would 
have to be made. Hovaver, we have 'Found that by making computations in the 
flow case at low Mach numbers we have been able to duplicate the no flow results . 
As an example, we have considered a linearly tapered, softwall transition section 
between two uniform hardwalled ducts. We have considered a wall admittance of 
A = 0.413 + i 0.720, a taper ratio of b g /b 0 = 1-268, a reduced frequency 
kb^ - 1.5, SL = 1.0 and b^ = 1.0. In Table 1 results are compared using the no 
flow program and the flow program with NF and PF modes. Using the 
PE modes we were forced to use M = 0,05 to avoid severe numerical problems which 
arise when the general equations are run at very low Mach number. This comes up 
because the leading terms in the v^ differential equations tend to vanish. 

This didn't occur using NF modes, but we feel that it is at least partly due 
to the difference in precision we used on the IBM 360 for the NF modes and the 
Burroughs B6718 we used for the PF inodes. In addition, for expedience , we 
used 80 integration steps and 3 modes in the PF case and used l 00 steps and 5 
modes in the NF case. Near zero Mach number the stability and accuracy of the 
integration is sensitive to the number of steps. Even with this difference in 


34 . 


Mach number we see that the results are definitely tending toward the no- flow 
results. These results support the accuracy of the terms in the general equations 
which are not flow related. They do not support a particular choice of basis 
functions because at small Mach numbers the NF and PF modes are nearly the same. 


(b) Convergence to the Hardwall case 

The computation of the transmission in a hardwall uniform duct provides a 
means to determine the accuracy of the uniform flow terms in the equations. To 
this end we have considered a high Mach number, M — -0.8, a nearly hardwall, 

A - 0.0001 + i 0.0001, at kb = 1.0, with b Q = b £ = 1.0 and l = 1.0. Both NF 
and PF modes were used. Table 2 is a comparison of results obtained analytically 
and results obtained from the knowledge that in a uniform duct no reflection 
or spurious mode generation occurs. The diagonal transmission terms are obtained 
by noting that for right moving waves the pressure solution is given by 

r„ -- A k e ‘ <y) 


The ratio of pressures at x = £ and x = 0 is 
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It is seen that these calculated results agree almost exactly with the results 
obtained computationally . Once again, this test case offers no insight into 
the choice of modes since in the hardwall case the NF and PF modes are the 
same. 

(c) Convergence to the Uniform Softwall Duct Case 

Transmission in a uniform softwall duct is also easily obtained by hand 
calculation and provides a convenient check case. Due to the particular for- 
mulation of our program we are only able to do this for the PF basis functions. 
We have considered the case M=-0.5, kb = 1.0, A = 0.720 + i 0.420, 
b Q = b^ = 1.0, Z - 1.0. The calculated results, obtained as in subsection (b) , 
except with defined from the computer program, are as follows: 

= /.^-237- £ / 0/03 A/,, -- 0.0/7/- £ 

-b 

k - 0.6787- i 4. ^ 008 S £ O 0067 
*2 

We have used 3 modes In the computational scheme and in this formulation 
essentially bypass the matching procedure which becomes unnecessary. 

Table 3 shows the results of the comparison and indicates complete agreement. 

Ihis check case does not verify the completeness of the PF modes simply 
because correct results are obtained. The PF modes are individually exact 
solutions to this problem. We are not asking them to produce a solution to a 
new problem. 

(d) Convergence to a Softwall, Uniform Segment Between Two Hardwall Infinite 
Ducts 

A somewhat more challenging problem from a computational point of view 
arises when we consider a uniform softwall segment between two infinite hard- 
wall ducts . The solution in the lined segment is exactly the same as obtained 
in the problem discussed in Sub-section (c) . However, in this case we are 
matching the solution to hardwall ducts rather than infinite softwall ducts. 
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Since reflection and spurious mode generation can occur we can observe tha 
performance of the NF and PF inodes, that is, we can observe if the solution 
in the lined section provides enough detail to predict the more complex 
acoustic field. To this end we have made two comparisons, one at low Mach number 
and one at high Mach number. 

In the low Mach number case we have used M = -0,1, A = 0,72 + i 0.42, 
kb = 1.0, b = b = 1.0, £ = 1.0. In table 4 we have compared transmission 

V A, 

and reflection coefficients using NF and PF modes with results from a mode 
matching scheme due to J.F. Unruh of the Boeing Co., Seattle, Washington. 

Because of the differences in the implementation of the mode matching scheme, 
only two reflection and transmission coefficients are directly comparable to 
our results. It is noted that the results of all three methods are nearly the 
same. At low Mach numbers the NF modes and the PF modes are not greatly 
different {they are identical at no-flow conditions) , so it is not totally 
surprising that they yield similar results. It is significant that the 
* results compare well with the mode matching approach, which is completely 

independent. 

As a second check case we have used M — -0.50, A = 0.720 + i 0.42, 
kb = 1.0, bg=b^=1.0, 1= l.O. We have again compared the current 
program using NF and PF modes with modtj matching results. The results of 
this comparison are given in Table 5. We notice a definite degradation in 
the comparison between the solution using the NF modes and the mode 
matching although the trends compare well. As we did not generate the mode 
matching results we can only assume their correctness at this point. 

Accepting this assumption, we feel that the degradation is due to the fact 
that the NF modes are not close solutions to the problem which involves flow, 
and hence it probably requires an increased number of basis functions to 
generate accurate results as compared to the no-flow case. The results 
using the PF basis functions are seen to be totally off. In Subsection (c) 
we used these basis functions to generate results for exactly the same 
conditions with the exception of the matching to uniform hard wall ducts 
. as in the present case. In Subsection (c) we obtained very good results, 

compared to exact hand calculations. Since the solution ir the "non-uni form 
, segment" is the same in both cases, it follows chat the poor correlation in 

the present case arises because of the matching. This indicates that the 
solution in the non-uniformity is not detailed enough to account for 
reflection and spurious mode generation. The original thought that the PF 
basis functions are not complete, and hence can't stand on their own as a 
■ — — ■ — ■— — — ■— 1 ■ — 1 ■ — J — ■■ — - ■■■ - — •«- • — - - — - j ~— — 
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aeries expansion for the softwall case, seems to be borne out. 

In order to more fully verify the trends when the NF modes are used, we 
need to develop an independent mode matching capability. The majority of 
the building blocks are available, but the press of other program development 
has prevented putting it together. 

(e) Comparison With a One Dimensional Formulation 

Davis and Johnson' ; have developed a one dimensional model of acoustic 
transmission through a compressible flow in a nonuniform hard walled duct. 

We have adapted this formulation to treat problems of the type encountered in 
the present investigation. The approach is to substitute an integration of 
their governing differential equation for the transfer matrix integration 
scheme in the program we have developed. The matching of the nonuniform 
segment to infinite uniform ducts is done in essentially the same way as 
explained previously. The Davis and Johnson (DJ) formulation can be considered 
a low frequency approximation, since only at low frequencies can we be 
assured of primarily plane wave propagation, as assumed in their analysis. 

To generate comparisons with the general program we have considered inlet 
(flow opposite to propagation) and exhaust (flow and propagation directions 
coincide) cases for a duct with a cosine shaped converging t t sr. In the DJ 
formulation the walls are hard and in the general formulation we take 
A = 0.0001 + i 0.0001. The reduced frequency is kb Q = 1.0, with = 1.0, 

= 0.75 and X. = 2.0. Mach numbers at x = 0 in the range -0.5<M^<0.5 are 
considered. These corresponded to Mach numbers at x = X in the range 
-0.93<M<0.93 . In the general program we have made the majority of our 
calculations using the NF modes but we have made runs using the PF basis 
functions at selected Mach numbers. Figures 2 and 3 show the comparison of 
results. The original NF results show good agreement with the DJ results 
for transmission coefficients for the entire Mach number range. Reflection 
coefficients show reasonable agreement in exhaust flows, but diverge from the 
one dimensional results for high inlet Mach numbers. The introduction of 
the PF modes has made a major improvement at high inlet Mach numbers and 
agreement is good, particularly for the reflection coefficient. We notiqe 
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some deviation in the transmission coefficient at M = -0.5, but in the 
light of the approximation in the method this is not considered sig- 
nificant. At this point the improvement is not thought to be due to the 
PF modes, which with hard walls should be the same as NF modes, but rather 
because of an improved implementation of the matching procedure at the 
ends of the nonuniformity. We consider this correlation to be of major 
significance in supporting the validity of the MWR. 


B. COMPUTATIONAL COMPARISONS USING COMPLETE FLOW BASIS FUNCTIONS 


To attempt to overcome certain shortcomings in the NF and PF basis 
function form of the Method of Weighted residuals, we have developed a form 
of the program utilizing the PNF basis functions. As noted previously, this 
form of the program is costly in terms of computer time and storage, but 
must be considered potentially more accurate because the basis functions are 
more nearly solutions to the problem, containing information on both positive 
and negative running modes. The PNF basis functions are probably complete, 
but they are not orthogonal, which creates coupling between the modes. 

In developing the PNF basis function formulation we encountered numerical 
problems which did not airse in our previous developments. In integrating 
the set of differential equations of Eqs. (37) to obtain the transfer matrix 
[T] , definied in Eq. (38), it was found that if the nonuniform duct segment 
is long (e.g. i = b Q ) , then ^^ents of the transfer matrix are large enough 
to make the subsequent matching operations numerically inaccurate. While 
this problem can be at least partially overcome by increased precision in 
the computations, it indicates a sensitivity that was not previously seen 
in equivalent cases in the original no flow problem or in the NF and PF 
formulations. We have circumvented the problem by considering only short 
segments (e.g. i = 0.1 b Q ) . This means that the analysis of a long duct 
will require the solution of a series of short duct segments and then a 
stepped matching procedure. In this situation we have essentially developed 
an acoustic finite element for a duct segment. We consider this to be an 
acceptable utilization of the method. We have also found the PNF case to 
be much more sensitive at low Mach numbers to the trend toward the break- 
down in the form of the differential equations caused by the duct mean flow 
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velocity approaching zero. 

A second major problem which has been noted is due to the use of both 
positive and negative running flow modes. The matrices [rfnml and L^nml 
whose elements are defined by 

-/ ^ ^7 ^y 

^ -- f<f>„ k d y 

o 

are apparently very nearly singular. We find numerically that the inverse of 
the leading matrix in Eq. (36) generates very large numbers (although our 
inverse scheme is known to produce an accurate inverse) which subsequently 
must be operated on through matrix multiplication on the right hand side to 
generate more moderately sized numbers. It is felt that this is a potential 
numerical trouble spot and the implications should be more carefully assessed. 
This appears to be telling us that even when the duct is soft walled and the 
flow does not vanish, the PNF basis functions are nearly linearly dependent. 

It may thus turn out that the right running modes and the left running inodes 
are not sufficiently independent to perform well computationally* We hope 
to be able to explore the implications of this observation in more detail. 

At the present time we are exercising caution in maintaining a double check 
on computational aspects such as the accuracy of inverses. 

In order to provide an indication of the correctness of the PNF program 
we have made several runs against cases which have a simple analytic solution. 
These check cases are: 

(a) Convergence to uniform hard wall duct results. 

(b) Convergence to uniform soft wall duct j ■r' -'ults. 

Table 6 shows the comparison of the hardwall duct results using PNF 
modes against the exact solution. We have used M = -0.5, kb = 1.0, 
b = b = 1.0, A = 0.0001 + i 0.0001, and % = 0.10. The analytical results 

are generated as in Subsection (b) of Section A. They yield 

¥ 

- 2.0 + t Or O TgA/J,. - o. 930/ - / 0./987 

K t 

U K ^ = O. &6G7- £ 3, 374 7~ F /l /7 22 - 0.7/2 0- £0.0475 

As shown in the table the computed and analytical results compare nearly 
exactly. 
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Table 7 is a similar comparison using results for a uniform softwall duct. 
In this case we have used M = -0.5, kb = 5.0, b = b =* 1.0, A = 0. 720+i0. 420, 

U A# 

and Z - 0.10. The analytical results yield 

7 u 

J<„ - 9. 7343- £ &&299S 7-&4AJ,, - <£>. ^ - / 0.&2 4 3 

/rS s 4. 9892 - i 6$ '3 7-<es?s/ 2 2 - 7/4Z - k & 6 O&S 

x 2. 

The comparison is seen to be very good. 

These comparisons, though limited in scope tend to support the implemen- 
tation of the PNF version of the Method of Weighted residuals. Because of 
the numerical problems that we have been forced to overcome we have not as 
yet mechanized the procedure by which we build up the solution for long ducts 
by mode matching rhe solutions for a series of shorter duct elements. This 
has prevented us from making the most demanding comparison against the Davis 
and Johnson low frequency, plane wave formulation. 

At the present time we feel that the PUP formulation is in basically 
good order. Check cases have been favorable and work is continuing toward 
implementing the mode matching technique. Sufficient experience has been 
gained to observe that the PNF form of the method of weighted residuals is 
computationally demanding and costly in terms of storage and run time. We 
are still working toward the point where the relative merits of PNF can be 
compared against the NF and PF implementations. 
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CONCLUSIONS AND EXTENSIONS 


The Method of Weighted Residuals in the form of a modified Galerkin method 
with weighted boundary residuals has been shown to be an accurate and efficient 
method for the analysis of the propagation of sound in nonuniform ducts with 
no flow. The method has been extended to include problems of propagation in 
ducts carrying a steady compressible mean flow. This extension must still 
be considered in the development stage in which extensive numerical experimen- 
tation is being carried out to optimize the implementation. Three different 
schemes are currently operational. These schemes differ in the type of basis 
functions used, the simplest implementation using modes from uniform ducts 
with no flow (NP) , the other two using positive running nodes, including the 
flow effect (PF) and positive and negative running flow modes {PNF) . 

The greatest amount of experience has been obtained with the NF basis 
functions. We have found that these modes yield good results over a wide 
Mach number range when compared against a simpler theory for low frequency 
propagation through a one dimensional flow in a hard wall, variable area duct. 

The method, as implemented, appears to break down at high inlet Mach numbers 
for calculation of the reflection coefficient. Results for transmission 
coefficient hold up through the entire subsonic range. Other check cases 
involving convergence to certain limiting cases of lined ducts have been 
generally good at low Mach numbers with degradation as Mach number increases 
(based on the presumption of accurate baseline results from other sources) . 

This is not unexpected since for ducts with soft walls the NF modes become 
increasingly poor approximations to actual solutions as the Mach number increases, 
and we expect to require a large number of basis functions to sythesiee a 
solution. 

The results with the PF basis functions have been mixed. In the hard wall, 
variable area duct our PF implementation gives good results over the entire 
Mach number range. This represents an improvement over the NF implementation. 
Since for hard walls the NF and PF modes ought to coincide, we currently 
feel that it is the implementation of the NF modes which the PF program has 
improved on. The improvement is most likely in the ...atching procedure. 

In soft wall ducts the PF modes are apparently not adequate. They fail in 'the 
limiting problem of a straight, lined section between two infinite hardwall 
ducts. It is felt that lack of completeness may be the problem. These modes 
are an exact solution in the lined segment. However, the transfer matrix 
generated when they are used is not detailed enough to account for reflected 
and spurious modes which arise due to the impedance discontinuity. Numerical 
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experimentation with PF modes is still in progress. 

The development of a computational scheme employing PNF modes introduces 
new computational problems not previously encountered. The problems have been 
identified and are due to the rapid growth in dimensionality caused by the 
doubling of the number of modes and the coupling between modes due to the lack 
of orthogonality of the modes. At the present time the technique for overcoming 
these difficulties is to consider only short duct segments, on the the order of 
1/10 of the duct height. This restriction does not render the MWR impractical, 
since longer nonuniformites can be treated by breaking the duct down into 
short sub-sections which can then be put together using mode matching tech- 
niques. When used in this context each sub-section can be viewed as an acoustic 
finite element. Because of the extensive development and analysis required in 
the PNF case we have only a limited computational experience. The comparison 
of results using PNF modes with certain limiting cases have been favorable. 
However, as noted elsewhere, computational costs rise rapidly because of the 
increased dimensionality. At the present time this implementation is still 
being developed and will be continued. 

We have accumulated sufficient experience with the three approaches to develop 
a strong preference for the use of the NF basis functions. This preference 
is based on the absence of severe numerical difficulties, the relative sim- 
plicity of the implementation and the relative expense when compared to the 
other two possibilities. It is believed that the problem noted in comparing 
to the Davis-Johnson one dimensional results can be eliminated so that the NF 
program will work as well as the PF program for this type of problem. However, 
we still recognize certain inadequacies. The most notable one is the apparent 
trend away from correlation with mode matching results as Mach number increases 
in the test case of a sc ft wall segment between infinite hard wall ducts. 

Our future research will be directed toward an improved form of the NF program. 
This will include a detailed analysis of the computations to identify problem 
areas. We currently favour a modification of the basis functions which may 
lead to an improvement in performance. It has been observed that the principal 
effect of flow on the basis functions is to change the lower modes, but leave 
the higher modes fairly close to the no- flow modes. H* seems possible that a 
set of basis functions using the lowest two or three PNF modes and the higher 
NF modes will improve the potential accuracy of the NF method. Of course we 
can't make an assessment of the completeness of this choice of basis functions 
without experimentation. 


As part of our future development we plan to write a simple mode matching 
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program to generate results as a baseline for comparison in our straight duct 
convergence studies. We have previously used results generated by other 
investigators, but we would prefer to develop independent results. 

Finally, in both the no flow and flow programs we have noticed a trend 
toward numerical difficulties if the basis functions include modes deep into 
cut-off. This problem is most severe when we are dealing with low reduced 
frequencies, since in this case it may require a number of cut-off modes to 
provide enough resolution to obtain a converged solution by the Galerkin method. 
We feel that it is the exponential character of the cut-off solutions that 
creates the difficulty. To alleviate this we currently have under development 
a modification of the computational scheme which essentially explicitly isolates 
this exponential behaviour so that the transfer matrix is not required to 
include it. This modification is nearly complete and if successful will be 
reported in the open literature. 
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FIGURE I 

DUCT GEOMETRY FOR EXAMPLE CASE 


COSINE TAPER ~ ty/b 0 =* 0.75 




Figure 2 


COSINE TAPER ~ bj;/b 0 = 0.75 



Figure 3 


Table 1 


CONVERGENCE OF FLOW CASE TO NO FLOW CASE 

bb_ =1*5 Linear Taper b./b_ = 1.263 A = 0.413 + 1 0.72 

0 ^ x 0 

% = 1.0 b Q = 1.0 

M=0.01 - NF Modes M=0.05 - PF Modes No Flow 

REFLECTION COEFFICIENTS REFLECTION COEFFICIENTS REFLECTION COEFFICIENTS 

12 12 12 

1-0.0776+10.2539 0.1431+10.0004 1 -0.0742+i0.2593 0.1504+10.0036 1 -0.0785+10.2524 0.1420-10 .0008 

2 0.0040+10.1522 0 .0340-i0, 0563 2 0.0175+10.1476 0.0303-10.0584 2 0.0009+0.1544 0.0343-10.0563 


TRANSMISSION COEFFICIENTS TRANSMISSION COEFFICIENTS TRANSMISSION COEFFICIENTS 

12 12 1 2 

1 -0.1628-10. 6314 0.0965-10.0073 1 -0. 1148-i0 . 6711 0.0973+10.0069 1-0,1741-10.6197 0.0960-10.0111 

2 0.1586+10.1949 0.1165-10.0326 2 0.1372+10.2082 0.1221-10.0309 2 0. 16 44+i 0.1910 0.1151-10.0332 




Table 2 


CONVERGENCE TO STRAIGHT HARD WALL DUCT 
kb = 1.0 b a /b 0 =1.0 A = O.OOOl+iO.OOOl, b 0 =1.0 l = 2.0 



:i i M=-0. 8 - NF Modes 

Selection coefficients 

| 1 2 

0 . 0000+iQ . 0000 0 . 0000-10 . 0000 

-0. 0003+i0.0002 0. 0001-i0.0001 

TRANSMISSION COEFFICIENTS 

:■ “0.8357+10.5446 -0. 0003+i0.0000 

[: -0 . 0001-i0 . 0004 -0. 0000+10 . 0001 


M=-0. 8 - PF Modes 
REFLECTION COEFFICIENTS 

1 2 

1 0 . OOOO+iO . 0000 O.OOOO+iO.OOOO 

2 O.OOOO+iO.OOOO O.OOOO+iO.OOOO 

TRANSMISSION COEFFICIENTS 

1 -0 . 8357+iO . 5445 O.OOOO+iO.OOOO 

2 O.OOOO+iO.OOOO O.OOOO+iO.OOOO 


n=-0.8 - Calculated 
REFLECTION COEFFICIENTS 

1 2 

1 O.OOOO+iO.OOOO 0. OOOO+iO. 0 A 00 

2 O.OOOO+iO.OOOO O.OOOO+iO.OOOO 


TRANSMISSION COEFFICIENTS 

1 ~0.8391+i0.5440 O.OOOO+iO.OOOO 

2 O.OOOO+iO.OOOO -O.OOOO+iO.OOOO 


Table 3 

CONVERGENCE TO UNIFORM STRAIGHT DUCT 
kb - 1. Q b^/b o = 1.0 A “ 0.72+10.42 b Q = 1.0 & = 1.0 


M~-0.50 £F Mode s 
REFLECTION COEFFICIENTS 



1 

2 

1 

0.0+10.0 

0 . O+iO . 0 

2 

0.0+10.0 

0.0+10.0 


M--0.50 - Calculated 
REFLECTION COEFFICIENTS 



1; 

2 

1 

Q.O+iO.O 

0.0+10.0 

2 

0.0' 10.0 

0.0+10.0 


TRANSMISSION COEFFICIENTS 


TRANSMISSION COEFFICIENTS 


1 2 


1 2 


1 0.0171-10.3635 0.0003-10.0001 

2 0. 0000+10. 0000 0.0081-10.0069 


1 

0.0171-10.3636 

O.O+iO.O 

2 

0.0+10.0 

0.0083-10.0067 



Table 4 


STRAIGHT 1 

SOFT WALL 

SEGMENT BETWEEN HARDWALL 

INFINITE DUCTS 

kb = 1.0 

b i /b 0 

A = 0.72+10.42 b Q - 

1.0 

Z ~ 1.0 

M=-0.1 " NF Modes 


M=-0. 1 - PF Modes 


M=-0 , 1 - Mode Matching 

REFLECTION COEFFICIENTS 


REFLECTION COEFFICIENTS 


REFLECTION COEFFICIENTS 

' 1 


1 


1 

I -0.195+10.167 


1 -0.1967+10.167 


1 -0.1922 +i0„ 170 

2 -0.030 +10. OBO 


2 -0.0287+10. 0806 


2 -0.028 +10.079 

TRANSMISSION COEFFICIENTS 


TRANSMISSION COEFFICIENTS 

TRANSMISSION COEFFICIENTS 

1 . 


1 


1 

1 0.187 “10.588 


1 0.1873-10.5882 


1 0.190-10.587 

2 0.057+10.042 


2 0.05712+10.0433 


2 0.059+10.0421 


Table 5 


STRAIGHT SOFTWARE SEGMENT BETWEEN HARDWARE INFINITE DUCTS 


kb = 1.0 

b/to =1.0 A = 0. 72+i0. 42 

b Q = 1.0 % = 1.0 

M=-0 . 5 - NF Modes 

M=-0 . 5 PF Modes 

M=-0.5 - Mode Matching 

reflection coefficients 

REFLECTION '.OEFFICIENTS 

REFLECTION COEFFICIENTS 

1 

1 

1 

1 -0. 10744-10.2030 

2 “0.0631+10.2742 

1 - 0.0798+10.3583 

2 0 . 0104+10 .1480 

1 -0.077+10.1665 

2 -0.091+10.2126 

transmission coefficients 

TRANSMISSION COEFFICIENTS 

TRANSMISSION COEFFICIENTS 

1 

1 

1 

1 0 . 0220-10 . 3353 

2 0. 0696-10. 0751 

1 0.0083-10.1211 

2 0.1590-10.2104 

1 -0.0103-10.3646 

2 0.057-10.084 



Table 6 


CONVERGENCE TO UNIFORM HARD WALL DUCT 
kb - 1.0 b^/b Q =1.0 A = O.OOOl+iO.OOOl b Q =1.0 £ = 0.1 

M=-0.5 - PNF Modes M=-0.5 - Calculated 

REFLECTION COEFFICIENTS REFLECTION COEFFICIENTS 



1 

2 


1 

2 

1 

0.0002+10.0000 

0. OOOO+iO. 0000 

1 

O.O+iO.O 

O.O+iO.O 

2 

0 . OOOO+iO . 0000 

0. OOOO+iO. 0000 

2 

O.O+iO.O 

O.O+iO.O 


TRANSMISSION COEFFICIENTS TRANSMISSION COEFFICIENTS 



1 

2 


I 

2 

1 

0,9798-10.1988 

0. OOOO+iO. 0000 

1 

0.9801-10.1987 

O.O+iO.O 

2 

0. OOOO+iO. 0000 

0.7120-i0. 04755 

2 

0 . O+iO . 0 

0.7120-10.04753 



Table 7 


CONVERGENCE TO STRAIGHT SOFT WALL DUCT 
kb = 5.0 =1*0 A = 0.720+10.420 b Q = 1.0 Z = 0.10 

M=-0.5 - PNF Modes M=-0. 5 - Calculated 

REFLECTION COEFFICIENTS REFLECTION COEFFICIENTS' 


1 

2 


1 


O.OOOO+iO.OOOO 

0.0002+10.0000 

1 

O.O+iO.O 

O.O+iO.O 

0.0000+10.0000 

0.0000-10.0003 

2 

0 . O+iO . 0 

O.O+iO.O 


TRANSMISSION COEFFICIENTS 

1 2 

1 0.5608-10.8243 0.0002+10.0002 

2 0.0000+10.0000 0.7144-10.6004 


TRANSMISSION COEFFICIENTS 

1 2 

1 0 ,5608-iQ . 8243 0.0+10.0 

2 O.O+iO.O 0.7142-10.6003 
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j A METHOD OF WEIGHTED RESIDUALS FOR THE 

j INVESTIGATION OF SOUND TRANSMISSION IN NON-UNIFORM 
| DUCTS WITHOUT FLOW 

■ W. Eversman and E. L. Cook 

j IVichita State University, Wichita, Kansas 67208. US.A. 

j AND 

J " R. J. BECKE.MF.YER 

' The Boeing Company, IVichita, Kansas 67210, U.S.A. 

I { Received 30 May 1 974. ana in revised form 18 July 1974) 

The transmission of sound in n non-uniform two dimensional duct without (low is 
investigated by a method of weighted residuals which leads to a set of coupled “generalized 
telegraphists* equations'*. Results for several duct configurations are compared with those 
from, respectively, a variational method, a stepped duct approximation, and an eigen- 
! function expansion method based on linearly tapered duct segments. 


! 1. INTRODUCTION 

I Relatively rigorous methods have been developed for the analysis and design of acoustically 
j lined and unlined ducts of uniform rectangular, circular, and annular cross section, with 

■ and without flow, and including the boundary layer effect in the flow case. The progress in 
I this area can be seen by referring to a few papers of the extensive literature which exists (see, 

for example, references The current capability in the mathematical modeling of 

duct propagation is limited primarily by the assumption of a uniform, infinite duct. 

There have been a Few recent studies directed toward the non-uniform duct probtem. In 
the case of ducts without flow a generally useful approach is the one developed by Zorumskt 
and Clark 16] for ducts of uniform area with lining variations and subsequently implemented 
1 by Alfrcdson [7] for the study of hard-walled ducts with varying cross section. This method 
j consists of representing the duct by a series of stepped ducts of uniform cross-section and 
j systematically accounting for the reflection and transmission process which occurs at the 
; intersection of the stepped elements. This procedure appears to be very useful provided it is 

■ used with due caution in the segmenting process. In the case of ducts with uniform area but 
. varying lining properties, it has been shown by Bahar [8]. for the case of electromagnetic 
' waveguides, to converge to the method developed in the present paper when the elemental 
! segments become vanishingly short. The principal difficulty with the stepped duct approach 
. as originally conceived is the high dimensionality of the numerical problem which results 

: A somewhat different formulation of the problem by the third author of the present paper 
; has reduced this difficulty while retaining the flexibility of the method, 
j Another method of general utility for ducts without flow is the variational approach of 
! Beckemeycr and Evcrsman (9). In this technique the acoustic problem is represented by a 
• variational principle and a Ritz minimization is employed to determine the coefficients in a 
, trial solution. The trial solution is in terms of basis functions which do not necessarily satisfy 
the boundary conditions (the boundary conditions are pan of the funci ional in the variational 
formulation) and do not have to be generated for each duct geometry. The stepped duct 
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‘ approach requires a set of eigenfunctions and eigenvalues in each stepped section and in the 
case of lined ducts this can be a problem in terms of computational requirements (the method 
introduced in the present paper must have eigenvalues and eigenfunctions at stations along 
the duct, but, as will be described, a rapid numerical scheme has heen devised for their 
computation). The variational approach sulfersa dimensionality problem in that complicated 
acoustic fields (both axial and transverse) require a large number of basis functions in the 
* trial solution. 

Other recent approaches to the problem are more approximate in nature by virtue of 
restricting geometry or frequency range allowed. Nayfeh and his co-workers have published 
several studies of propagation in non-uniform ducts with and without flow. The paper by 
Nayfeh, Telionis and Lekoudis (10] is representative. They restrict themselves to ducts with 
; slowly varying cross-section, lining properties and flow properties and employ a perturbation 
' scheme. To within the level of accuracy which they retain, they do not predict the generation, 

; reflection, and transmission of modes other than the one incident on the non-uniformity. 

, It appears that a higher level of approximation is required to predict this. 

Karamcheti and his co-workers have also made contributions in this area. King and 
: Karamcheti [II] studied plane waves in ducts with one dimensional flow by a method of 
, characteristics and Huerre and Karamcheti [12] used a short wave approximation for the 
: same type of problem. Similar problems were studied earlier by Powell [13] and Eisenberg 
and Kao [14]. 

Tam [ 1 5] seems to have published tlte first paper dealing with a multi-modal approach to the 
problem of non-uniform ducts with flow. His technique is a perturbation scheme based on 
the assumption of slowly varying cross-sectional area. The first order approximate solution 
is obtained by Fourier transformation. 

A recent paper by Cummings [1 6] studies the novel problem of the acoustics of a wine 
bottle. The wine bottle is a non-uniform duct without flow and the acoustic field is approxi- 
' mated by a Rungc-Kutta integration scheme based on the Webster Horn Equation. This 

■ method allows only a plane wave mode of propagation and hence is limited to the lower 

■ frequency ranges. There have been many studies of the horn equation since Rayleigh first 
' introduced it and it is a favorite topic in texts on acoustics. 

In the present research program we are interested in the multi-modal propagation ofsound 
in non-uniform ducts of fairly general shape. The finul goal of the program is the study of 
propagation in non-uniform ducts with flow, so that we will be interested in methods for the 
no flow case which appear to be extendable to the ease with flow. Of the two generally applic- 
j able techniques mentioned previously, the variational method does not seem to be readily 
. extendable to the ease with flow. The stepped duct approximation might have some applica- 
; tion, although it is certainly questionable whether the non-uniform flow field can be repre- 
sented in sufficient detail in a scries of stepped uniform segments. We are thus led to look for 
' another method with promise for the flow case. In this paper we will introduce the method 
and assess its utility in the ease without flow, as there are equivalent results available against 
which a direct comparison con be made. In the course of the development and application 
of the method it has become apparent that the method of weighted residuals is an important 
alternate method for the duct with no flow, and indeed, may well be superior to the other 
’ two methods or general utility. 

The method or weighted residuals (MWR) employed here actually was first employed in 
connection with electromagnetic waveguide problems by SchclkunofT [17, 18), The field 
equations for these problems arc identical to the classical acoustic equations in certain cases. 
Schclkunoff's work followed work by Stevenson ( 19, 20] which used MSR but led to a some- 
what diTere'nt formulation which does not seem to have been widely used, Stevenson appears 
. to be the first to suggest the application of methods of this type to acoustic horn problems. 
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Bahar and his co-workers studying ionospheric propagation of microwaves as a terrestrial 
waveguide problem have made extensive use of MWR (see, for example, reference [8]} and 
Reiter [21] has formalized the approach. 

j In this paper we have used MWR to approach the problem of multi-modal propagation 
in non-uniform ducts without flow. We deal with a two dimensional duct in which both area 
variations and lining variations are permitted. The results obtained by using MWR arc com- 
pared with those, respectively, from a stepped duct theory, the variational approach, and a 
j segmented duct theory in which the duct segments arc radial (sectoral). In reviewing the 
; literature it does not appear that results for lined ducts of variable cross-section in the multi- 
! modal case, with fairly general area variations permitted, have previously been presented. 

2. METHOD 

; In this analysis two dimensional ducts of infinite length arc considered. The extension to 
circular, annular, or rectangular ducts is straightforward, but of course more demanding 
. computationally. The extension to include finite duct terminations is easily included in the 
; present formulation provided reflection and transmission characteristics of the termination 
| are available. 



Figure 1. Duel configuration. 


Figure 1 shows the type of configuration under consideration, in which two semi-infinite 
, uniform duct sections of wall impedance z, and z c and areas s t and s c arc joined by a transition 
section of length / of variable area ^(.-e) and variable impedance z»(x), where x is the axial 
. co-ordinate. In this analysis the area variation will be restricted to be continuous, but the 
lining variation can be discontinuous at the ends of the non-uniformity. 

In terms of dimensional acoustic pressure, p \ panicle velocity, V*, and density, p*. the 
acoustic equations far complex harmonic motion of the form e ,a " in time are 

; kup* + podiv V* ~ 0, 

iiupo V* ■» -grad p*, 

; p* - c 5 p*. 

. where the ambient state is given by p 0 , p 0 , V 0 = 0 and c « (ypolPoY 11 is the ambient speed of 
. sound. By introducing the non-dimensional variables p “ p^JpoC 1 , p «* p*/po, V *• V*/c and 
eliminating p one can write 

ikp + divV**0 
itV — — gradp 

where k = wjc ** 2 nfX and it is the free space wave number. 


(I) 

( 2 > 


.4 



In order to specify a boundary condition at the duct walls, a co-ordinate system as shown in 
J Figure 2 is introduced. This figure shows the manner in which the duct height profile and 
1 slope of the height profile are specified as well as the designation of the local outward unit 
! normal at the duct wall, ». The duct wall boundary condition employed is characteristic of a 
normally reacting lining in the presence of a harmonic pressure variation: namely, 


j P *-2 K;-zV*-v, (3) 

where Yf is the component of the acoustic particle velocity in the direction of the outward 
normal, z is the wall impedance which may bea function of axialposition. In non-dimensional 
' j variables equation (3) becomes 

t s ' 

I «?«“— V-v. 

PaC 


V-v- Ap, 


! at the duct wall, where A - p a cfz is the acoustic admittance ratio of the lining. In the case of 
a uniform duct of two dimensional, circular, annular, or rectangular cross section equations 
t (1), (I) and (4) can be combined to produce the classical problem or propagation in a lined 
! duct which has been thoroughly stud icd. In these problems powerful techniques can be brought 
: to bear since a co-ordinate geometry can be chosen to make both the field equations and the 
. boundary conditions variables separable. For general non-uniform ducts this is not possible 
t and an alternate approach must be used. 



Figure J. Duct geometry for example problems. 


To demonstrate the method of weighted residuals (see Fintayson*s book [22] for a com- 
plete description of the general method), we consider a two dimensional duct as shown in 
Figure 3 in which the wall al y ■=> 0 is hard and the wall at y »* is lined with a material of 
admittance /If*). This representation is also valid for a duct symmetric about the x-axis 
with the acoustic propagation also symmetric with this axis. The governing equations are 
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equations (1) and (2) with the boundary condition obtained by expanding equation (4): 

v ** 0, y - 0, (5) 

pcos0“irsin0K/fp, (6) 


j where a and v arc the non-dimensional axial and transverse velocity components and tan 9 *■ 
< dbldx is the wall slope. 

| 7/c seek solutions to this problem in the form of a finite scries of specified basis functions: 


1 

Ph" f P*(x) cos K M y, 

m-i 

CO. 


1 u.(x)coSK m y, 
■»1 

(8) 

1 

2 tV(*)sinxvy, . 

W 


J where the. * « are defined as the infinite sequence of eigenvalues or 
j Kbtmxb — McbA. 

] The basis functions are recognized as the eigenfunctions for propagation in a uniform duct 
, which has the same height and admittance as the non-uniform duct has locally. This means 
| that the eigenvalues k. are functions of x and that the basis functions change with axial 
position. Note that These basis functions do not individually satisfy the boundary condition 
of equation (6); however MWR will force them to do so collectively. The choice of basis 
1 functions may seem to be unnecessarily complicated since wc will be required ro provide a 
1 new set of eigenvalues and eigenfunctions at each axial position; however, we will show that 
' this can be done quickly and easily. Any disadvantage is outweighed by the more rapid 
i convergence of this type of technique when the basis functions are chosen to represent as 
1 nearly as is practical the actual solution. 

' If the trial solutions, equations (7H^>. are substituted in the differential equations and 
I boundary conditions, they do not in general satisfy them, but instead leave an error, or 
' residual. Thus, in terms of these trial solutions one can write 



1 /£* «<p*cosf>— u*sin0- Ap*),^. 

Since a continuous function must be zero if it is orthogonal to every member of a complete 
set we will force these residuals to be orthogonal to every member of the complete set formed 
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by the basis functions themselves at every station x along the duct axis: 
- * 


j ^i(x.y) cos.K.ydy = 0, 

D 

00) 

> 

J -Ri(x.y) cos x M y dy - 0, 

O 

0D 

» 

J* iJj(x,y) sin K.ydy - 0, 

O 

(I2)_ 

Jt((x,£}cosic # 6 — 0. 

(13) 


i The operations indicated by equations (10)-(13) yield, for *** 1,2,. ..AT, 

' Ut f p*COSK,ydy + I — ^cosK,yd^+ — -£O5K«ydy«0, (14) 

I i { dx l dy 

r f dPn 

ik J u w cosK«ydy+ J — costf.ydy **0, (IS) 

' r f dpN 

; ik J % sin + J ™sin ^ m ydy - 0, (16) 

(u^cosO- u w sin(J- -4 /j„)cosk,6“0. (17) 

I • 

By using Leibnitz's rule for differentiation of integrals containing a parameter the partial 
derivatives of u H and p M with respect to xin equations (14) and (1 5) can be replaced by ordinary 
J derivatives. For example: 

: f^jcosK.ydy ~~ Jtr w cosK„ydy- J ^ [cos k, y] dy - u w(x, 6) cos k, A — 

1 o o o 

The partial derivatives with respect to yean be eliminated by integration by parts: 

/ • »g p t 

—cos ic^y dy - o^x, b) cos k, b + tc. J v„ sin x t y dy. 
dy 5 

O 

| * 

sinK.ydy-op^x.&jsin^b-ff, fpNCOSK.ydy. ( 18 ) 

i i dy . * 

Equation (16) can be used to eliminate r. from equation (14) and the boundary residual, 

, equation (17), is used to simplify the boundary terms in the resulting equation. In this way, . 
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1 equations (14) through (17) can be written in terms of two equations: 

j 5x/ Uh 005 K * ydy ** ~ lfc [ 1 ~ (■£") j Jp « 005 K *> d y + / u «“(cos *r,y) dy - 

[ i*&4lp*(x,b)cosK„& 

Tb (l9) 

j p N cosK Jl ydy-*-ii|M„cosK»ydy*hJp w ~(cosK,y)dy-fp A <x,6)coSK,6^. (20) 


j The trial solutions for p N and ti H are substituted into equations (19) and (20). In the case 
■ of the duct without flow the eigenfunctions are orthogonal so that 


m 

J cosK.ycosic^ydy- 


1 where 6 m =*Q,£&m, and &„ = 1. 
Equations (19) and (20) become 


I where 


I""*! 1 ( J t)] P ‘ , ~ „?i “ Ji ULP "' 

I PLjS,^U2...Jf. 
dx «-i 

. . h»» 

1 dAL I dv, r . 

u: -“wd7^ + FdrJ 


i kbA 

COS X, fc COS K m b 

cost? 

N m kb 


tan 0 

J Ft, “ 1C, - -jj— cos K m bcosK m b. 

I Equations (21) and (22) can be represented in matrix form by 



— L 

‘HtI 

<5- 

-WM 

w 




-if; 

L) 

. 

r. 


, It is noted that if the duct is uniform so that dx^dx^ Oand di/dx- tan 0«*O then equations 
i (21) and (22) reduce to equations for the axial field in a uniform duct corresponding to each 
1 mode of propagation. In this case they arc a form of the telegraph ists' equations. When the 
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duct is non-uniform, coupling terms arc introduced and the resulting equations arc “generaj- 
ired telegraphists' equations", a term introduced by SchclkunofT (17], 

liquations (21) and (22). or their equivalent matrix form in equation (23), will define the 
acoustic field in the non-uniformity providing one can specify initial and or final conditions 
on the//, and p. at the start of the non-uniformity at jr = Oand the endx — /.This can be done 
by assuming that the sound is incident from the left, or x < 0, and the acoustic field for * < 0 
1 consists of incident and reflected waves. For x> l, in this analysis, we assume the sound 
propagates to the right in a semi-tnfmilc uniform duct, and hence consists only of waves 
. propagating away from the non-uniformity. In general, one can write, for x<0, 

i p{x, y) - 2 fat c *' 1 * 1 + on c'H cos R.y, 

m 

where al, an are coefficients of the incident and reflected waves, respectively. arc the eigen- 

values for propagation in the uniform duct for x < 0, and 

T • 

where the plus sign is chosen if kjk is real and if fcjk is complex the sign is chosen to make the 
imaginary part negative. One can write similar expressions for the uniform duct, x > l: 

/ 

a(x,y) - e “ u: * “ 4je , *S*]cosi£j/. 

p(x, >-) - 2 [bt * + bn e 11 * 1 ) cos y, 

where the starred quantities apply for x> i. At jt = 0 one can match particle velocity and" 
pressure: 

- /£] cos if. 2 u,casK M y, 

£ fat + aT] cos y - 2 p„ cos K,p, 

U 1 the same number of uniform duct eigenfunctions as there arc basis functions in the non- 
uniformity is used, then the orthogonality of the eigenfunctions can be used to obtain a matrix 
relation between the oj and an and the u, and p.at x = 0. This is 


■*»] [( tH-#) 8 " t , 


where, at x = 0 
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r 

, cos K.y cos x m ydy “ iV. 5^ 


H 


cosK t ycosit m ydy. 


An exactly analogous relationship can be obtained at x ** /: 


pJJ) PL \ PL 


! c *»r t j ft* 


r where, atx*/, 


? 

, — J COS K, 


m ycosK„ydy~N K S m 


Z* 

PZ M = j cosK.yooiKZydy, 

Hence, it is possible to obtain relations between the wave amplitudes in the uniform ducts 
and the pressures and particle velocities in the non-uniform section in the form 


H--R 


| By using an integration scheme such as a fourth order Runge-Kutta, equation (23) can 
; be integrated from * «=* 0 to x «* / to obtain a transfer matrix relating u„ p M at x - / to u t and 
5 p m at a;t=*0: 


; By using equations (25) and (26), the wave amplitudes for x > I can be related to those for 
ar<0: 


n- 
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£ 77 ^*] - [TA n ] - ITA^lTAnY'lTAnl ( 31 ) 

The physical significance of the elements of these matrices is as follows: 

REF„ •» Amplitude of reflected mode f due to incident mode j with amplitude 
1-0+ iOO. 

TRAN lt “ Amplitude of transmitted mode / due to incident mode / with amplitude 
1-O + iOO. 

In the formulation in the preceding equations we have referenced both reflected and trans- 
miticd waves with respect to x <= 0. Because oflargc exponential factors which can occur for 
cut-off modes in the e terms in equation (24), it appears to be more appropriate to reference 
the transmission coefficients to x = /. In this case, one can rewrite equation (24): 



All of the steps leading to equations (30) and (31) can be repeated and a new set of reflection 
coefficients and transmission coefficients can be derived. The reflection coefficients will not 


ill 


v 
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' change but the transmission coefficients referenced tore/ will be related to those at x = 0 

: by 

1 [TRANl. 0 -[ e Z'][TIUNl„,. 

I Once the reflection coefficients are obtained, one can write the initial values for uJLx) and 
I P*(x): 



If the pressures and particle velocities in the non-uniformity arc required, equation (23) can 
be integrated with equation (32) being used as initial conditions. In our compulations wc 
have stopped with the reflection and transmission coefficients for the non-uniformity, 
principally due to the extra computational time implied to obtain the complete acoustic field. 

One can also easily account for termination conditions for x > l other than the semi- 
infinite duct employed in thi* analysis. If a relationship is known between bi and b7 at the 
end of a finite uniform duct, i.e., 

5 {*>-}“ !*]{&*). 

’ where [RJ is the termination reflection matrix, then equation (29) can be written as 



Then 


[ ' isnbn-ir^.ifa^+rr^iKo-} 

) and 

i {c-} — (t-RHT^] - - {TXi.lHo*}. 

! m - IlTVfii) - [7Vi n K[*1l7Vf»l - [7Vl„r‘([K]i™ u ) - f7Vl a , ])]{<»*}• 


3. IMPLEMENTATION OF THE METHOD 

j The most serious potential disadvantage of this form of MWR is the need to compute the 

, eigenvalues which appear in the basis functions and hence in the coefficients in equation (23), 
j As noted earlier the eigenvalues, k„ arc roots of 


1 


k/j tan ni *« i kb A. 


(33) 


i The usual technique for obtaining the eigenvalues is a Ncwton-R.iphson iteration based on 
| suitable initial guesses. Wc expect to require 5 to 10 basis functions to obtain solutions 
sufficiently converged for our purposes, and hence wc wilt require 5 to 10 eigenvalues for 


' each point required in the integration scheme used to generate the iraurir; matrix. The 
J iteration scheme is likely to be too costly computationally to be of practical value. Wc have 
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circumvented the problem by replacing equation (33) by the dificrent’al equation which it 
satisfies: 


1 + cos(2sfc) 

4x 2nb + sin (2Kb) 


/I dir 1 dd\ 

2i kb A I-— + ] . 

\b dx A ebej 


We start w ith values of Kb at x = 0, for each basis function required, and generate the required 
values of Kb at each integration step of equation (23) using a Runge-Kutta integration of 
equation (34). Our scheme is configured to start with no initial guess other than hard-wall 
eigenvalues foreach required basis function atx = 0. This approach hasproven to be far more 
rapid than iteration and has shown good accuracy as can be checked by occasionally evaluat- 
ing Equation (32) and noting the size of the residuals which develop. 

The coefficients which appear in equation (23) have previously been given. The specific 
terms in these coefficients arc given below; 

+ sin2*„f>\ 

•“n 2 Zb )• 

dN„ I d6 sm2tf.fc — 2n m bcos2K.bdK m b 
dx bdx " Qxtb) 1 dx 

/ * . , lfsitt («» — **) 4 sin(ic,+ 0^1 

b f cos (k\ ~K m )b COS (K, + K m ) b\ 

2[ (*„ - K m ) (a'. + O J 


m 

jy sin K t y cos K m y dy ■ 


6 1 fs:n2K B 6 cos 2 k. b 


2 (2*,h) 1 2k. b * 


[ , , tAMlcoSK lt i>co5»r„f> ( I \cosic.bcosK m b a 

k. b tan k. b -j — — I -- I n 

cosflj N m kb (cost) J N. 

These expressions, together with equation (33) ore used to generate U&,, and P^ M , 


4. RESULTS 

1 Very few results have been published for the multi-modal transmission of sound in non- 
uniform ducts. The work of Zorumski and Clark [6], a sequel by Lansing and Zorumski [23], 
i and the work of Alfrcdson [7] arc the most nearly equivalent to the present method. They 
used the stepped duct approach for segmented linings or area changes in ducts with an open 
end. Since we have concentrated on the non-uniformity and have chosen not to represent 
terminations, we have no common basis for comparison. The approximate theories mentioned 
previously are considered too limited in scope to make comparison feasible when the labor 
and expense of doing so is considered. The only extensive results which we have been able to 
’ use conveniently are those of Beckcmeycr and Eversman [9] who developed a variational 
. approach and compared it with stepped duct approximations and an approximation in which 
the duct was segmented into linearly tapered sections. Of course, that study wasclosely related 
to the present one and the thrust of the entire research program has been to develop a basis 
of comparison. 
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J As a first example we consider linearly tapered hard-wall transition sections joining two 
| uniform ducts. Forthesc cases a 15° taper transition is used. The first example is a diverging 
| taper from a duct of heigh! b 0 to one of height (I -Man 15")fr 0 ,thelcngthofthenon-uniformity j 

! being equal to its initial height. The second example is a converging taper w hich iscxactly the 
’ reverse of the diverging one. We have compared these results with those generated in reference [ 

j [9) and reproduced here as Figures 4-7. In the diverging ease kb 0 values of t-Q, 2 0, 2-5, 3-0 i 

j and 3-5 are compared and in the converging case kb, values of 1*0, 1-5, 2-0, 2*5, 3-0 arc used. '■ 

At kb — 2-5 in the small duct the second mode is just cutting on in the large duct while at 
kb = 3*5 in the small duct it is cut on throughout. The results from reference [9] for the diverg- t 

ing taper include those from two levels of stepped duct approximation (5 sections and 10 
sections having been used, respectively), from the variat ional approach, and from a segmented 
duct theory (eigenfunction expansion) in which a linearly tapered duct segment with solutions 
given in terms of Bessel functions was used. In the last method pressures and velocities were 
matched at the ends of the taper, x «=* 0 and x = /, by collocation. No account was taken of the 
slight difference between the planar co-ordinate surface at die ends of the uniform duct and 
the circular co-ordinate surfaces of the radial segment. Details of the characteristics of these l 

methods can be found in reference [9]. In the converging taper only the 5 section stepped 1 

duct and the variational approach were used. The M WR was used with 6 basis functions and ' 

50 integration steps axially, ( 



it> 0 


1 Figure 4. (a) Transmiuion coefficient and (b) reflection coefficient for mode 2 with mode I incident. 
Diverging 35” linearly tapered unlined trantition section, bt ** I -27 fro IJ"; L — A a. — — 10 Sccuon 

stepped duct; , S section stepped duct; a, Fourier bases; a, eigenfunction expansion; □, weighted 

residuals. 

| Figures 4{a) and 5(b) show the transmission and rdleclion coefficients for the diverging 
| duct Tor the second mode with the first mode incident at x — 0 with amplitude 1*0 -i-O-Oi. 
I We have made the most detailed comparisons for the second mode since with increasing 
j reduced frequency it goes through the cut on phenomenon and hence is considered the most 
i challenging computationally. The transmission coefficients, shown in Figure 4<a), are in 
1 good agreement throughout the frequency range, all methods considered yielding essentially 
the same result. In Figure 4(b), it is seen that the results for the reflection coefficients for 
j kb 0 > 2*5 arc not in as good agreement. However, it is noted that the variational method, the 
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Figures. (a)Tranimisiion coefficients and (b) reflection coefficients for modes I and 2 with mode I incident. 
Converging IS* linearly tapered unlined transition section. 60 = l-27fti.a=* 15°;2.«6/. 


Mode 1 

Mode 2 

Method 

( 



5 section stepped duct 

A 

* 

Variational (Fourier bases) 

□ 

O 

Weighted residuals 


eigcnfuiiciion expansion method and the MWR remain in close agreement but they do not 
agree closely with either stepped duct result at the higher kb 0 . The imaginary part of the 
reflection coefficient is the first to show significant deviation. We have concluded that our 
choice of the number of modes in the stepped duct approximation was not adequate at the 
higher frequency, particularly after the second mode begins to cut on. 

| Figures 5(a) and (b) combine the results for transmission and reflection coefficients in the 
first and second modes with mode ( incident for the converging taper. In this case the 5 section 
stepped duct and the variational method are used Tor comparison. Agreement is good through- 
out the frequency range. The stepped duct approximation appears to be more adequate in 
this case than in the diverging case, particularly for reflection coefficients, since more modes 
are used on the reflection side than were used in the diverging duct case. 

, As a second comparison we have analyzed al seleclc; frequencies the transmission and 
reflection characteristics of a lined diverging transition section between uniform hard-wall 
ducts. Wc have used a cosine shaped transition which has the same small and large heights 
and the same length as the linear taper previously discussed. The only available results for 
comparison arc from the variational approach. We encountered convergence problems with 
the variational scheme as our limited computational capacity permitted neither double 

' precision nor more tha" 40 input basis funciions. A scries of results were obtained for which 
4 axial cosine waves. 4 axial sine waves and from 5 to 9 transverse functions were used. In all 
cases convergence was observed to be occurring with increasing modes, but in several cases it 
was not reached to our satisfaction (judged by comparingsuccessivc runs with more and more 


modes). For this reason wc have only made comparisons w here convergence in the variational 
' case was acceptable. The MWR was used with 6 transverse modes and 50 axial integration 
■ steps. Wc have generated results in the form of reflection and transmission coefficients in 
modes 1 through 3 for mode 1 incident. Results are compared in Table I. Absence of varia- 
I tional results indicates a lack of convergence at the level of approximation used. Wherever 
| converged results are available, the MWR shows acceptable agreement with the variational 
! approach. Because of the computer limitations noted in the variational results, the MWR 
| results are considered the more accurate. 


j Table I 

i Reflection and transmission coefficients for a cosine shaped lined transition section, 
1 First mode incident; amplitude » 1 -0 + 0-0i 


kb a 

Method of weighted residuals 


Variational method 

r 

Reflection in 
first three nude; 

Transmission in 
first three modes 

i 

Reflection In 
first three modes 

Transmission in 
first three modes 

1*5 

-0 05 4- 0*27i 

— 0*18 — 0-Gli 

—003 4* 0 2Si 

—0 17 - 0 6J* 


-0 00 + 0-14i 

017 + 0 20i 

001 + 0I3i 

0 I 8 + 0 20i 


0 01 -0 031 

— 0 04 — 0 03i 


— 0 04 — 0 04i 

20 

0-15 + 0 03i 

-0-26— 0-40i 

0*17+ 0 06i 

-0 25— 04 li 


0*I5 + 0 21i 

0*54 + 0 02i 

016 + 0l7i 

0 55 — OOOi 


-0 01 - D-05i 

-0 0B + 0-03i 


-0 10 + D02i 

2-5 

— 0-11 — 0 02i 

—0-47 — 0 44i 

-O-ll -0 021 

-0 48 - 0 44i 


0- 19-0 021 

0 23- I 02i 


0 19- 1-021 


-0 02 — 0-04i 

Ool +0 08i 

-002-0 031 

0 01 + 0 lOi 

30 

-0-03 + 0 04! 

— C 67 — 0*141 

-O 03 + 0O4i 

—0 67 — 0-1 3i 


0 25-OOli 

-0-41 — 0-58i 

0 I9-005i 

— 0 42 — 0 55i 


— 0 02 - 0 04i 

002 + 0 05i 

-0-02 - 003i 

002 + 006i 


; Mote: at each reduced frequency the first two columns are the conversed mutt from MWR. The second 
two columns arc from the variational method when a converged result is available. 

kb a =- 1‘5, A “0-413 + 0 720i kb t •« 2 5. A « 07SS + 0* J36i 
kb a ~ 2 0, A = 0-664 + 0'720i kb a - 3 0, A - 0760 + 0 MOi 

i The results with which we have compared ore all at relatively low frequency al which at 
* most two modes propagate. Wc do not currently have available higher frequency results 
■ against which to compare, principally due to limitations in the computer facilities which we 
| have used. However, we can gain a measure or confidence in the MWR by investigating its 
I convergence properties at higher frequencies, We have considered the cosine shaped diverging 
transition section at a reduced frequency of kb 0 = 12-57 and fora lining with A ■= 0 76 + 0 30i. 
1 In this case in the small duct the plane wave and 3 higher modes arc well cut on and mode 5 is 
, just cutting on. In the large duct the plane wave and 5 higher modes propagate (kb, =• 1 5*94). 

We have used four combinations of basis functions and integration steps beginning with 6 
' transverse modes and 50 steps and progressing through 8 modes and 50 steps. 6 modes and 
| 100 steps, and 10 modes and 50 steps. Increasing the number of integration steps made 
virtually no difference. Increasing the number of basis functions provided a very definite 
' convergence trend. To illustrate this convergence of the method wc have shown in Table 2 
i the reflection coefficients in mode 5 due to all incident modes (REF U ). Mode 5 was chosen 
'. because the driving reduced frequency is right at the mode 5 cut-on frequency in the source- 
side duct. This causes relatively large reflected components in this mode for any incident 
: mode. In addition, in Table 2 wc have shown the direct transmission coefficients in all incident 
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Table 2 

Convergence trends in cosine shaped diverging transition. 
kh n . 12-57. A - 0-76 + 0-31i, b 0 = 1-0, b K = 1-26 8.L=1'0 

Reflection coefficients in mode 5 Direct transmission coefficients 

Incident a , , * * 

mode 6 jj^j, p -^5 g Bas ; s pnig 10 Basis FNS 6 Basis FNS 8 Basis FNS 10 Basis FNS 


1 

-0223 

-0215 

-0 2122 

0-736 

0736 

0-736 


+0 IS2i 

+0 I41i 

+0-!38i 

+0034i 

+0 034i 

+0034i 

2 

0230 

0-222 

0-2186 

0-523 

0523 

0 522 


— 01 Sli 

— 0 I39i 

-0 136i 

+0 213i 

+0-21 2i 

+0-21 2i 

3 

-0-245 

-0-232 

-0-228 

0-146 

0-148 

0-148 


+0-1 46i 

+0 I32i 

+0-1 29i 

-+0-45H 

+0-4511 

+045H 

4 

0-260 

0-238 

0232 

-0-337 

-0-341 

—0342 


— 0-191 i 

—0-1901 

— 0-I88i 

+O-024i 

+00291 

+0-03H 

5 

—0 900 

-0896 

-0 896 

0016 

0017 

0018 


-LO-OSOi 

+0 0801 

+0 080i 

+0-012i 

+0012i 

+U012i 

6 

0-279 

0 272 

0269 

-0-057 

-0 060 

-0 061 


-HH521 

+0-1611 

+0 I64i 

— 0-012i 

-0016i 

-00171 

7 


-0194 

-0 191 


-0-005 

-0 005 



-0-1 28i 

— 0 I30i 


+0-0021 

+0 002i 

8 


0-157 

0-154 


-0-002 

-0 002 



+0-108i 

+0 11 li 


+00021 

+0 001i 

9 



-0-130 



-0 001 




— 0-097i 



+0 00H 

10 



0-114 



-0-000 




+O087i 



+0-0011 


; modes {TRAN a ). Coefficients arc shown for 6, 8 and 10 basis functions. The sequence of 
results is distinctly convergent. In this case the use of 8 basis functions is probably sjfficicnt 
from an engineering standpoint. 

As a final example of the use of the method we have computed the acoustic power balance 
for the linear converging taper between uniform hardwall ducts in the case when the taper is 


t Table 3 

I A caustic power balance in linear taper converging transition section. 


kb 0 - 1 5-94, b 0 = 1 -268, 6, « I -0. L => 1 0 


Incident 

mode 

Incident power 

Reflected power 

Transmitted power 

Absorbed power 

A •* 0-76 + 0-30i 

l 1 000 

0 003 

0754 

0-242 

2 

1 000 

0 007 

0595 

0 399 

3 

1000 

0010 

0 672 

0-317 

4 

1-000 

0015 

0-571 

0414 

5 

1000 

0 020 

0-307 

0674 

6 

1 000 

0 169 

0050 

0-781 

A “ 0-0001 + OOOOli 
1 1 000 

0004 

0996 

0-000 

T 

1-000 

0 0075 

09923 

0-000 

3 

1 000 

0-026 

0 974 

0000 

4 

1-000 

0-148 

0852 

0 000 

5 . 

1 000 

0557 

0443 

0 000 

6 

1-000 

0948 

0052 

0-000 
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lined (A = 0 76 + 0-30i) and also in the hard-wall case at kb 0 = 1 5-94. For each propagating 
mode incident on the non-uniformity we have computed the incident power, the reflected 
power, and the transmitted power carried in all resulting propagating modes. The power 
dissipated in the acoustic lining is the amount by which the sum of the reflected and t ransmitted 
powers fail to match the incident power. Table 3 shows the results of these computations for 
both cases. The hard-wall case (A = 0 0001 -fO OOOIi) is significant in that it adequately 
approximates the known result that no power is dissipated in the non-uniformity. While wc 
have not at this point made extensive parametric variations, it appears from these results that 
the reflective characteristics of a non-uniformity can be enhanced for modes near the cut-off 
condition. We intend to report more detailed investigations of specific non-uniformities at a ' 
later date. 


; * 5. DISCUSSION AND CONCLUSIONS 

The method of weighted residuals has been shown-to be En accurate and rapidly convergent 
method for the compulation of the acoustic transmission and reflection properties of non- 
uniformities in duct; without flow. Other methods applicable to the problem have been 
implemented and numerous comparison runs have been made to validate the MWR and to 
build confidence in its use. In all cases where other results could be generated, comparable 
results were obtained. As with any numerical technique, care must be exercised in the use of 
MWR. The choice of basis functions appears to favor good convergence characteristics, but 
in any new physical situation convergence experimentation should be undertaken. 

One of the reviewers of the original version of this paper offered the suggestion that basis 
functions which satisfy 


\p' + k* d' «=* 0, tji'ffi) > 




be considered. This is equivalent to using an effective admittance A * »• A/cosO. If this is done 
the eigenvalue equation is 

! . , i*M .... 

Kb tan Kb . (35) 


The array of coupling coefficients, I/'*, then vanishes. While this is attractive it should also 
be noted that this would complicate the computation ofdtcjdx. and hence i/J.and P ’ m . and 
effectively introduce more coupling there. Hence, whether this is really a better set of basis 
functions would require numerical experimentation. The important point to be made is that 
the optimum use of methods of this type depends on extension experience in the particular 
application. 

It has also been pointed out that the differential equation for the eigenvalues, equation (34). 
becomes singular when the impedance assumes the Cremcr optimum value (2]. This is 
because the optimum is defined by a double eigenvalue and, as Zorumski {24] points out, a 
double eigenvalue implies that not only is 


K m btanx.b - ikbA 


! but also 


■ [ir, b tan K m b - \kbA] ■» 0. 
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This means that 

2k, b + sin 2 a, b = 0 

at a double point. When this occurs it is apparent that equation (34) is singular. Of course 
this problem is not unique to this method. The Ncwton-Raphson iteration would display a 
similar problem in computing eigenvalues as a function £>r some changing parameter if the 
double point were encountered. If this occurs and if it cannot be circumvented by testing for 
the occurrence and employing an expansion of the differential equation at the singular point 
to gel past it, then a different set of eigenfunctions might be employed. For example, those 
generated from equation (35) would not have the same double points as ones generated front 
equation (33). 

Our philosophy has been that with improved analysis methods for finite, non-uniform ducts 
with multi-modal propagation, the practical importance of the Crcmer optimum has been 
reduced in that the optimum lining will generally not correspond exactly to the double point. 
We view the potential difficulty as important but not crucial to the utility of the method. 

The non-uniformities treated in this paper arc of the type which one might encounter in 
applications. They are not small, nor are they very abrupt {although all of the results for 
lined ducts presented here have a c scontinuous change from hard-wall to soft-wall), n is to 
be expected that the more radical the non-uniformity, the less rapid will be the convergence, 
and the mo rc basis functions and i ntegration intervals will be required to achieve a sal isfactory 
result. This is an inherent property of this type of method. At this point we arc satisfied that 
we can treat problems or practical importance. 

The results to date conclusively show that it is important to treat the duct non-uniformity 
problem from a multi-modal standpoint. A given incident mode will generate spurious 
reflected and transmuted modes which can have an important effect on acoustic lining design 
and radiation properties. In addition there appears to be a possibility of using geometry 
changes and their attendant reflective properties to enhance acoustic attenuation, but the 
multi-modal performance of the duct will have to be known for design purposes. 

The problem of the most immediate importance is the extension of the method to the case 
of ducts with flow. In contrast to other methods investigated, the MSR appears to be extend- 
able w hen flow is ; ■'solved. Of course the complexity of the problem expands considerably, 
but the basic numerical scheme can be modified and expanded to include this case. Our 
research program is currently directed toward achieving a workable method when flow is 
present. 

After the original version of this paper was written an expanded version of the work of 
Zorumski and Clark [6) has appeared as a NASA document [25]. The workinthisdocumcnt 
is quite approp: .ate to the development of the interface relations of equations (24a) and (24b) 
and is closely related to the development c r eem ; ms (25) and (26). 
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COMPUTATION OF AXIAL AND TRANSVERSE WAVE NUMBERS FOR 
UNIFORM TWO-DIMENSIONAL DUCTS WITH FLOW USING A NUMERICAL 
INTEGRATION SCHEME 

1. INTRODUCTION 

; The purpose of this letter is to detail a scheme by which the axial and transverse wave numbers 
‘ for propagation of sound in uniform two-dimensional ducts with uniform flow can be com- 
puted. In the method used the eigenvalue equation is transformed into a first order non- 
! linear ordinary differential equation. By using appropriate initial values this differential 
I equation is integrated by using a Ruogc-Kutta algorithm to yield solutions which are the 
transverse wave numbers for the duct. The transverse wave numbers then arc used to compute 
1 theaxial wave numbers. The method proposed is particularly useful for the rapid computation 
of a number of duct eigenvalues at a single reduced frequency, lining admittance and duct 
flow Mach number. It has the advantage of being an inherently stable computational scheme. 
In addition, the ordering of the eigenvalues is well defined by their relationship to a cor- 
responding set of hard-wall eigenvalues (the initial values in the integration scheme). 


2. method 

The specific problem is that of propagation of sound in a uniform two-dimensional duct of 
height b with one wall hard and the other with a lining admittance A. This configuration can 
be viewed as modelling a duct with this lining arrangement or a duct of height 26, sym- 
metrically lined, with symmetric propagation. Prop.'.gatio* at reduced frequency kb is 
considered, where k is the plane wave wave number k * w'c, to being the driving frequency 
and c the speed of sound. The duct Mach number, M, is assumed to be subsonic. 

It is well known that for this situation propagation in the duct can be represented by a 
1 superposition of acoustic modes of the form 

p m - A. c'“' c~ w '-* cos k. y, 

where x is the duct axial co-ordinate and y is the transverse co-ordinate measured from the 
hard wall. k„ is the transverse wave number /or the mh mode of propagation and k Mm is the 
corresponding axial wave number. They are defined by the simultaneous eigenvalue equations 



i AH 


f 

il -M 



(D 




For the present analysis equations (1) and (2) arc considered in the somewhat modified 
for ns 


kb tan £6^ j - 
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1 + M 1 1 -(1 - A/ 2 ) 


1 — A/ 3 ' w 

Throughout the analysis it is important to maintain consistency In the ± sign choice. The 
significance of the sign choice is straightforward for ducts with attenuation. In this case the 
axial wave number is written as 

kjk » a ± if}. 


r , /«vt° 
- (l -"'>(*,] ■ 


and the square root is the principal value taken in the upper half plane (positive imaginary 
pan). The modes of propagation are in the form 

/>, = /f»e“"'"*> r ‘ , e t/ * 3 ‘cosR._y. 

In order to have attenuation with increasing x it is necessary to choose the negative sign in 
equation (2) or equation (4). Acoustic modes defined according to this convention will 
correspond lo propagation in the positive x direction. If the flow is taken to be in the positive 
x direction, tbv. this will represent propagation with the flow (exhaust mode). When the 
, positive sign is chosen in equation (2) or equation (4) one defines propagation in the negative 
x direction, against the flow (inlet mode), In the case when p is identically zero, which occurs 
when xjk =0 or njk is real and 1 - (1 - A/ 2 )(k 7A) 2 > 0, the principal value of the square 
root is taken on the relative real axis. (This is done so that the sign choices given above always 
relate the same way to right and left running waves.) 

The usual way to obtain selected eigenvalues of the doubly infinite sequence of and k„ t 
defined by equations ( I ) and (2) or equations (3) and (4) consisis of using an iteration scheme 
based on suitable initial guesses. Users of this method arc well aware of potential instabilities, 
difficulties encountered in ordering the results, and the requirement for accurate initial guesses. 
Successful numerical schemes based on the iteration approach arc in existence but generally 
have to be quite sophisticated to overcome the problems indicated. The method that intro- 
duced here is both simple and accurate. 

One considers the eigenvalue tejk to be a function of some parameter, ij. In general, one 
also would consider k b, A and A f to be functionsofthis parameter. If equations (3) and(4)are 
differentiated with respect to tj and combined, then the following single ordinary differential 
equationresujtjS: 




d .. _ 2iAw 


-MS 


and r* 2 is the principal value or the square root as defined previously. 

If n'k ts known for some combination A, kb, M. then equation (5) is a differential equation 
defining the variation of k k as A , kb, m change :t% functions of ij from the original values. 



Equation (5) can be used in two different ways. In the simplest application, and the one which 
is the major subject of this paper, one considers kb and M to be independent of ij. Then the 
differential equation describes the variation of Klk as A varies with rj. Hence, if A is given some 
prescribed variation, the corresponding variation in x/k can be obtained by integration. 

Consider the problem of finding the eigenvalues defined by equations (3) and (4) for a 
specific condition A, kb, M. If A = 0 the eigenvalues are known to be 


k 


(«-!)« 
kb * 


n •= 1, 2. 


( 6 ) 


Now let A vary according to A(t}) - tjA f , where 0 < tj < I. Equation (5) becomes 



If equation (7) is integrated on 0 < tj < 1 , a herd-wall eigenvalue from thcscquencc ofequation 
(6) being used as the initial value, then at tj = I the solution to equation (7) will be an eigen- 
value for ihe condition A f , kb. A/. Furthermore, it can be identified as the eigenvalue which 
branches from the hard*wall eigenvalue on the real axis In the x]k plane which was used as 
an initial value. Each hard-wall eigenvalue on the real axis branches into two soft-wall 
eigenvalues, one for propagation in the positive x direction and the other for propagation in 
the negative x direction (corresponding to the sign choice), Hence, starting with .V hardwal! 
eigenvalues, one can compute 2 N eigenvalues (n exhaust mode eigenvalues and .V inlet mode 
eigenvalues) for the duct conditions A f , kb, A f simply by integration of equation (7) on 
0<ij< 1 with hard-wall eigenvalues as initial values. For each eigenvalue, xk, one can 
compute the corresponding axial wave number, k,lk, by using equation (2), 

A second application of equation (5) has proved useful in the study of non-uniform ducts; 
in this case A, kb, and M can be considered to be functions of the axial co-ordinate, x. In this 
case, equation (5) defines the variation of the eigenvalues, xlk, with respect to axial position. 
The starting values in this case would be the eigenvalues corresponding to the conditions at 
the starting point, say x= 0. The eigenvalues at this starring point could be computed con- 
veniently by using equation (7). 

Yet a third application, not implemented as yet by the author, would provide a means of 
carrying out parametric variations to determine the variation in duct attenuation sws with 
independent variations in A, kb or A/. For example, if it were desired to compute the variation 
in attenuation for a given lining admittance at a given Mach number for variations in kb 
one could construct a simple functional dependence for kb, say 

kb{n) « n ikb),, 

{kb\ t being the highest kb required. If the initial values are chosen as know n eigenvalues when 
i ; takes on its initial value, rj„, and if the integration is carried out on jj 0 < t) c I, then one 
rapidly can generate the results of the parametric variation. 


3. IMPLEMENTATION | 

The computational scheme follows nearly exactly the theoretical development described f 

above. A fourth order Rungc-Kutia integration package has been used. In solving the basic t 

eigenvalue problem as defined by equation (7), a variable step size has been used, starting 
with 2y„ of the full interval for the first five steps and then switching to 10!. for the last nine 
steps. To possibly increase the accuracy of the results a Ncwton-Raphson interaction has 


1 - 
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been used at the end of the integral ion for each eigenvalue to refine the result which may have 
some accumulated error due to the rather coarse integration grid. Experience to date shows 
very little need to do this. The best way to judge the error is to insert the computed eigenvalue 
at the current integration step into the eigenvalue equation, equation (3), and evaluate the 
residual developed. It has been found that if a significant residual develops it usually will be 
in the first step away from the Initial value. It may prove worthwhile to use the developing 
residual to automatically vary the step size. 

The only instance in which equation (7) or equation (5) requires special attention is when 
it becomes singular. This always occurs in using equation (7) when one S.'-ips away from the 
hard-wall eigenvalue, K{k = 0. It also can arise when a double eigenvalue occurs. While the 
secund situation is possible, no special precautions to account for it have been taken by the 
author, the philosophy being that the chance oHt occurring ts remote unless one is specifically 
attempting to compute the double eigenvalue. It is worth mentioning that the useofan auto- 
matically variable step size based on the current residual provides an excellent means of 
controlling the accuracy when the integration passes near a singular point. 

When it Is required to step away from the hard-wall eigenvalue, K}k = 0, in equation (7), 
the first step is made by noting that for small A and small k I k one can write 


K 

k 


\kbAAf 

T 


(I ± Af), 


where A is the initial step away from r/ = 0 which in the author’s implementation has been 2% 
of the total interval. This result is then refined, by using a Newton-Raphson scheme, and used 
as a starting value for an integration beginning at fj »= 002. 
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